最近着手把CSK移植到DSP中,先看一些DSP中图像处理的一些例子,第一件事当然就是怎么把图像数据倒入CCS工程中了,去年倒是用过一点CCS,再拿起来已经忘得差不多了,这篇文章主要记录一些学习的过程:
目录
- DSP导入图像数据
- 窗函数和加窗
- 定点数和浮点数
- 图像的FFT运算
- 图像的IFFT运算
一. DSP导入图像数据
搞了一下午大概可以了,主要是如何导入数据,如何利用CCS的Image analyzer来做显示。
1. 准备dat形式的数据
CCS导入数据是dat格式的,利用matlab把数据保存到dat文件里,这个直接找的代码:
clear
clc
I=imread('img.bmp'); //这里把图片都进来就行
[s1,s2]=size(I);
Num=s1*s2/4;
fid=fopen('img.dat','w'); //“写”形式打开
fprintf(fid,'1651 1 80000000 0 %X
',s1*s2/4); //信息头
for i=1:4:(s1*s2-3)
fprintf(fid,'0X% 02X% 02X% 02X% 02X
',double(I(i+3:-1:i)));
end
fclose(fid); //关闭
不用去纠结得到的dat的格式规则,要用到的话再细究。
2. 准备内存,导入数据
可以利用CCS的load memory功能,这个功能是在调试的时候可用的,所以先开辟内存,然后再load memory。这里也比较简单:
load memory
点击next选择内存位置和长度:
load memory
注意这里数组不能太长,可能会造成溢出,我把堆栈设置改成0xFFFF,这里用2500的长度就可以了,具体这块还没有搞清楚,50*50的图片也够我测试了。点击finish就可以完成导入了,然后可以在memory browser里查看数据:
memory browser
需要设置的参数我都用红 {MOD}标出了。
查看图像。
查看图像用tools->image analyzer这个工具,会弹出两个框,一个属性框properties设置一些参数:
参数设置
然后在对应的图像窗口点refresh就可以看到图像了:
image
加载图像加载进去倒是可以,但是是访问不了的,写了个很简单的阈值处理跑不了。
替代方案
最后也没有发现到底是什么样的原因,但是在李老师(李翔宇)的指导下,把图片数据的首地址给定,最后给在的是0x00850000,这样就可以了,老师说也 不知道什么原因,先放下吧,最起码这样可以进行一些简单的算法仿真了。
具体这样:
unsigned char *img=(unsigned char *)0x0850000;
再load memory就没有问题了,也可以进行操作。
二.窗函数的实现和导入
CSK的实现过程中要用到两种窗函数,分别是高斯和汉明,这两种窗函数可以利用matlab提前生成好,然后作为头文件来导入到CCS工程中。这个实现起来也不难。
matlab代码
%得到汉明窗和高斯窗的c代码,不能直接写入h文件,就先写入txt再复制过去了,主要是要中间的逗号。
fid=fopen('cos.txt','wt');
cos_window = hann(64) * hann(64)';
cos_window=cos_window(:);
fprintf(fid,'cos_window_64_64={');
for i=1:size(cos_window)
fprintf(fid,'%f,',cos_window(i));
end
fclose(fid);
这样生成的是txt文件,主要是包含逗号就可以了,稍加修改就可以复制到CCS建立的头文件中了,数据类型选择float类型。我试着导入了一个汉明窗,64*64的,拉成一维数组后导入是这个样子:
汉明窗
CCS中如果想看二维图像长什么样子的话还得把float转换成8位int,挺麻烦的,所以就这样看看吧。
然后我们把窗函数加在图像上看一看长什么样子,我这里换了一副图像,所有的操作都是针对64
64的图像来的。
导入的图像为
原图
加窗是要把两个数组相乘,直接用循环把两个数组乘起来就可以了(element-wise),floatuchar=float,这样得到的结果如下:
img_with_window
然后就是把乘的结果转换为uchar型的来显示,首先定义一个存放转换之后结果的数组,然后用循环逐一转化(这种应该都是可以用多核进行优化并行计算的)我一开始是这么写的:
for(i=0;i<4096;i++)
{
*(img_smooth++)=*(img_with_window++); //转换为uchar
}
会莫名奇妙出现问题,不仅img_smooth的结果不对,还会把原先img_with_window的数据给冲掉。数组地址我都是手动给定的,保证不会冲突的,也没找到什么原因。(如有人知道还望指点!)
img_smooth
img_smooth就变成这样了,这肯定是不对的。而且img_with_window就全变成0了。
后来换了一种方法用下标来做,就没有问题了,真是神奇。
image_smooth
转换为uchar之后长这样,对比可以知道这个是正确的,这个时候就可以看加窗之后的图像了:
image_smooth
都是没有问题的,这个smooth的命名的原因就是因为上次写均值滤波时用的那个变量,还没有改。
三.定点数和浮点数的区别
PC编程很少遇到这么细节的问题,但是DSP上就不同了,以前只知道定点数需要定标,浮点数是采用类似于科学计数法的一种方法,具体的细节就不清楚了,DSP还有定点和浮点之分,所以把这里的细节看了看,然后整理下。
定点数
所谓定点数就是约定小数点的位置是不变的,这样通常将定点数据表示成纯小数或者纯整数,纯小数的话,就把小数点固定在数值部分的最高位之前;纯小数的话,则把小数点固定在数值部分的最后面,如图:
小数点在机器中是不表示出来的,是一个提前约定好的位置,对于一台计算机,一旦确定了小数点的位置,就不再改变。
必须明确的一点就是计算机是无法计算小数的,所有的小数计算都是通整数计算完成的,这就导致事先约定小数点的位置尤为重要,这就是定标,在定点运算中,定标很重要。
为了简单,我们用16位的来说吧,也就是说参与运算的数都是16位的整型数,那么如何处理小数的关键就在于如何确定一个数据小数点的位置,这样的话数据的精度和数据的范围就是一对矛盾了,精度越高,范围则会越小。
数的标定有Q和S两种表示方法:
对于16位的来说,最高位为符号保留位,那么根据定标的不同,有Q0到Q15共16中定标方式,精度依次增高,范围依次降低。
例:
0010 0000 0000 0000b 表示8192 Q0定标
而同样的一个数:Q15定标
0010 0000 0000 0000b 表示0.25
这个图片不是特别清楚,可以看下:
DSPlib中的fft32x32函数就是以Q15定标的,所以如果要用这个函数,就要对数据进行转换。这个可以自己手动定标,也可以用DSPlib里的函数,我找到一个函数叫做DSP_fltoq15(),这个是float到Q15的一个转换函数,是可以进行这样的转换的。这个在CCSdsplib的文档中就有。
这个算法是直接给出的,也很简单,实际上是一个移位(乘法就相当于移位),注意调用之前x应该是归一化到[-1,1)之间的。
void fltoq15(float x[], short r[], short nx)
{
int i, a;
for(i = 0; i < nx; i++)
{
a = 32768 * x[i];
// saturate to 16.bit //
if (a>32767) a = 32767;
if (a<.32768) a = .32768;
r[i] = (short) a;
}
}
至于浮点数怎么判断溢出和运算规则就暂且不细究了,我这次的应用只涉及到所用的算法是要求Q15格式,所以需要手动转换。
浮点数
正是由于定点数的表示方法太过僵硬,固定的小数点位置决定了固定位数的整数部分和小数部分,不利于表示特别大和特别小的数,现在绝大多数的现代计算机系统都采纳了浮点数表达实数。
浮点数用一个尾数,一个基数,一个指数来表示一个数。
比如123.45用十进制的科学技术法可以表示为:1.2345*10^2。其中1.2345就是尾数,10是基数,2是指数,这样可以更加灵活得表示更大范围的实数。
IEEE浮点数
IEEE标准规定了两种基本的浮点格式:单精度和双精度,其中单精度具有23位有效数字,总共占用32位,那么自然是有8位指数,1位是保留符号位,双精度共64位,其中有52位尾数,11位指数,1位是保留符号位。
单精度
双精度
现代的X86架构早已摒弃定点数,在可以进行浮点数进行计算的计算环境中,可以自由得利用float表示小数进行计算。
细节不说了,一般不是很涉及底层的话也用不到这些知识。
四.图像的FFT运算
我在另一篇
文章里总结了,二维的FFT可以转换为一维的FFT进行计算,即先对每一列进行FFT运算,然后再第一次FFT运算的结果上进行列FFT变换,具体的证明去看那篇文章,这样二维的就可以用一维的来做,利用的是dsplib的fft32x32函数,这个函数可以快速计算FFT,先来看这个函数:
fft32x32
这里的w是提前生成的,利用dsolib里的一个exe文件就可以了。
exe
用cmd来打开这个:
直接看截图吧,这个是一位老师指导的,使用说明也给的很清楚,这样之后就会生成一个文件,这个文件可以是头文件,也可以是c文件,我是直接生成c文件,把系数复制到我的程序里了。
另外x[2
nx]和y[2nx]都是32位的数据,虚部紧跟实部连续存储,32位数据采用Q15定点,所以这里涉及到一个数据转换,整个fft二维算法的大框架我是借鉴我们组另外一个老师的,所以不便放出来,需要把数据转换成Q15格式的,这里的转换建议先归一化到[-1.1),然后左移15位赋值给int型的,这样就能得到Q15的数据。
nx是FFT变换的点数,无需多说。
看下结果:对于一个64*64的图像,我用matlab和ccs分别进行了计算,先归一化,然后去做FFT,这里只看实部:
matlab
ccs
精度是有损失的,但是看着前几十个数的结果还算对的,但是我把这个结果导入到matlab里和matlab计算的fft的结果相比(我是把数据拉长做了个差),发现一个问题:
实际上是每四行的第一行结果是正确的,其他行的结果是有问题的,这个我估计是它的频谱是乱序的,具体什么原因现在还不清楚,等我下午再去研究研究,应该是顺序出了问题。
我先把中间的结果输出,即每一列先做FFT的结果保存起来看一下,我发现这个结果是在正确的,在matlab里和上面同样做了差,两位小数精度下是完全准确的。
所以说这个fft32x32的函数是不要求输入乱序的,这个我翻手册的时候也没发现要求输入乱序。我试了试暂时还没有找到结果错误的原因。我也试了下看是不是乱序,以matlab_fft变换的第二列为标准,在ccs的结果中找最相似的一整列,实验证明并不存在这样的一整列。说明不是乱序的结果。
搞了半个下午,还是没有发现问题所在,我去干干别的再来搞这个,是在不行要去请教下老师!
18年1月19号上午更新:
刚才正在看别的东西,突然灵机一动,想了想是不是内存地址互相冲突呢?于是赶紧打开CCS试了一下,第一次没成功,赶紧检查发现还是有一个变量的地址没有完全指定,导致两个变量的地址是完全一样的,我的程序里是xtemp和ytemp,指定了以后赶紧去看第二行第一个数,果然就和以前的不同了,导入matlab简单画了个图,perfect!
matlab和CCS傅里叶变换
误差基本在两个三个小数点之后,可以接受!
五.图像的IFFT运算
FFT成功之后,还是花了半个下午的时间把测试的代码整理成函数放到头文件里,因为代码不多,就没有新建源文件,之所以花了半个下午是因为当时一个内存死活写不进去数,对于一个int型的数据取float直接就是零,搞得我都想把电脑砸了,后来重启了一下CCS竟然好了,神奇的错误,这也是现在不想做太多底层的原因,特别是和硬件相关的。
这两天爸妈来西安了,中午刚送走,从车站赶回来刚赶上下午上班时间,实际上都困得不行了,就想去睡觉,到了办公室把早上改的IFFT看了一下,其实IFFT就只需要把FFT的代码中的dsp_fft32x32改成dsp_ifft32x32就可以了,其他的东西都不用改,因为其他得辅助代码只是为了先行后列这样操作罢了。
但是算的结果确差了很多,早上没想到底是怎么回事,下午回来把数据导入到matlab里reshape显示一看,竟然是原图,那么肯定是数据整个扩大了一个整倍数,找了一组对应的数字一除,4095.6,简直完美,这才想起逆傅里叶变换是有个系数的,DSPlib中的ifft函数并没有包含这个系数,因为是行列做了两次,所以应该是1/64*64=1/4096,于是在最后除上这个数就ok了。
这样这个CSK算法最难的地方就过去了。
2018/1/24更新
六:C674Xdsplib导入及使用
昨天碰见老师了,问我做的怎么样了,交流了一下决定先用浮点的DSP来做,如果速度可以的话再移植到定点上去,所以前面用定点做的要改一改,于是又去下载了C674x的daplib的库,这里面有一个dsp的函数是浮点的,还有一些矩阵乘法数组乘法之类的东西,昨天一天都在鼓捣这个东西,因为网上的资料实在是太少,搞得有一阵子都想把电脑砸了,不过赶在晚饭之前还是弄好了,不过就是这样也很气,感觉做了许多无用功,没人带做东西真的累。一晚上做梦也都是浮点定点这点破事,一天好像都在看内存里存的什么东西。
不过抱怨归抱怨,还是把这里总结下,如果有人能看到且要用到,希望有用。
DSP官方出了许多库函数,前面用的dsp64x的库函数,打开之后解压之后长这样:
用的时候,我当时是直接把lib里的四个文件全部拷贝到项目中来,这里面长这样:
全部复制过去,然后项目里吧include文件夹import进去,用的时候要包含相应的头文件,这个现在不用了不说了,很简单。
重点说下dspC674x的dsplib的调用,这个方法可以用在许多lib的导入上,比如mathlib,DSPc674的dsplib需要用到mathlib,所以这两个都导入了,首先去下载两个库,TI的官网上就可以直接下,exe格式的,安装完成后是一个文件夹,长这样:
readme里说了这些东西都是干嘛的,主要的东西都在packages是里,这个里面是:
下面说配置:
CCS,项目的properity里:
把这几个都添加进去。
还有一个:
linker里面把lib添加进去,就可以用了,用的时候只需要添加#include"dsplib.h"就可以了,这个头文件打开看一下就明白了,它实际上是一个链接头文件,包含了dsplib里所有的函数,当然也可以用哪个包含哪个。
这样准备工作就完成了,这个浮点库里做FFT用的是这样一个函数:
点开可以看原图,输入输出都是float形式的,虚部紧跟实部保存,N是数目,这里面主要有两个参数需要解释,第一个是
brev是一个乱序表,第二个是ptr_w是旋转矩阵。
brev用这个就可以了,也有说这个参数根本就没用,我还是写上了。
unsigned char brev[64] = {
0x0, 0x20, 0x10, 0x30, 0x8, 0x28, 0x18, 0x38,
0x4, 0x24, 0x14, 0x34, 0xc, 0x2c, 0x1c, 0x3c,
0x2, 0x22, 0x12, 0x32, 0xa, 0x2a, 0x1a, 0x3a,
0x6, 0x26, 0x16, 0x36, 0xe, 0x2e, 0x1e, 0x3e,
0x1, 0x21, 0x11, 0x31, 0x9, 0x29, 0x19, 0x39,
0x5, 0x25, 0x15, 0x35, 0xd, 0x2d, 0x1d, 0x3d,
0x3, 0x23, 0x13, 0x33, 0xb, 0x2b, 0x1b, 0x3b,
0x7, 0x27, 0x17, 0x37, 0xf, 0x2f, 0x1f, 0x3f
};
另外有一个函数:
void tw_gen (float *w, int n)
{
int i, j, k;
double x_t, y_t, theta1, theta2, theta3;
const double PI = 3.141592654;
for (j = 1, k = 0; j <= n >> 2; j = j << 2)
{
for (i = 0; i < n >> 2; i += j)
{
theta1 = 2 * PI * i / n;
x_t = cos (theta1);
y_t = sin (theta1);
w[k] = (float) x_t;
w[k + 1] = (float) y_t;
theta2 = 4 * PI * i / n;
x_t = cos (theta2);
y_t = sin (theta2);
w[k + 2] = (float) x_t;
w[k + 3] = (float) y_t;
theta3 = 6 * PI * i / n;
x_t = cos (theta3);
y_t = sin (theta3);
w[k + 4] = (float) x_t;
w[k + 5] = (float) y_t;
k += 6;
}
}
}
我想直接调用这个函数出不来,觉得还不如自己生成数据拷贝进来方便,FFT和IFFT用的是一套w,所以只需要生成一次复制进自己的头文件或者源文件就行了。我是直接在VS里把这个函数复制进去,生成的数据拷贝到CCS的头文件中,竟然出奇得顺利就完成了。
未完待续---
2018.7.16: 这个暂时就终结了,这个算法最终大部分的函数都写完了但是没有调通,主要的原因可能还是在我比较抗拒DSP吧,所以做的成果都打包发给老师了,老师让一个职工往下做了,希望可以有些不错的结果。