DSP图像处理

最近着手把CSK移植到DSP中，先看一些DSP中图像处理的一些例子，第一件事当然就是怎么把图像数据倒入CCS工程中了，去年倒是用过一点CCS，再拿起来已经忘得差不多了，这篇文章主要记录一些学习的过程：

一. 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\n',s1*s2/4);   //信息头
for i=1:4:(s1*s2-3)
    fprintf(fid,'0X% 02X% 02X% 02X% 02X\n',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

需要设置的参数我都用红色标出了。

查看图像。

查看图像用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，挺麻烦的，所以就这样看看吧。

然后我们把窗函数加在图像上看一看长什么样子，我这里换了一副图像，所有的操作都是针对6464的图像来的。

导入的图像为

原图

_加窗是要把两个数组相乘，直接用循环把两个数组乘起来就可以了（element-wise）,float_uchar=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文件，把系数复制到我的程序里了。

另外x2_nx和y2_nx都是32位的数据，虚部紧跟实部连续存储，32位数据采用Q15定点，所以这里涉及到一个数据转换，整个fft二维算法的大框架我是借鉴我们组另外一个老师的，所以不便放出来，需要把数据转换成Q15格式的，这里的转换建议先归一化到[-1.1),然后左移15位赋值给int型的，这样就能得到Q15的数据。

nx是FFT变换的点数，无需多说。

看下结果：对于一个64*64的图像，我用matlab和ccs分别进行了计算，先归一化，然后去做FFT，这里只看实部：

matlab

ccs

精度是有损失的，但是看着前几十个数的结果还算对的，但是我把这个结果导入到matlab里和matlab计算的fft的结果相比（我是把数据拉长做了个差），发现一个问题：

实际上是每四行的第一行结果是正确的，其他行的结果是有问题的，这个我估计是它的频谱是乱序的，具体什么原因现在还不清楚，等我下午再去研究研究，应该是顺序出了问题。

我先把中间的结果输出，即每一列先做FFT的结果保存起来看一下，我发现这个结果是在正确的，在matlab里和上面同样做了差，两位小数精度下是完全准确的。undefined

所以说这个fft32x32的函数是不要求输入乱序的，这个我翻手册的时候也没发现要求输入乱序。我试了试暂时还没有找到结果错误的原因。我也试了下看是不是乱序，以matlab_fft变换的第二列为标准，在ccs的结果中找最相似的一整列，实验证明并不存在这样的一整列。说明不是乱序的结果。

搞了半个下午，还是没有发现问题所在，我去干干别的再来搞这个，是在不行要去请教下老师！

18年1月19号上午更新： 刚才正在看别的东西，突然灵机一动，想了想是不是内存地址互相冲突呢？于是赶紧打开CCS试了一下，第一次没成功，赶紧检查发现还是有一个变量的地址没有完全指定，导致两个变量的地址是完全一样的，我的程序里是xtemp和ytemp,指定了以后赶紧去看第二行第一个数，果然就和以前的不同了，导入matlab简单画了个图，perfect！undefined

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吧，所以做的成果都打包发给老师了，老师让一个职工往下做了，希望可以有些不错的结果。