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

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏每日一篇技术文章

OpenGL ES _ 高级03_调试工具

11310
来自专栏吾爱乐享

白盒测试的测试方法及基本路径测试法

25130
来自专栏LET

谈谈3D Tiles(2):数据结构

49240
来自专栏贾老师の博客

一致性哈希 Consistant Hash

14250
来自专栏YoungGy

MMD_3b_StreamAlgorithms

Stream概述 DBMS的区别 Stream模型 Query种类 应用 Sliding Windows 简介 例子 Bloom Filter motivati...

17890
来自专栏轮子工厂

教你用翻译软件快速阅读大量英文文献

对于一些引用的英文文献,我们需要快速地了解整篇文献讲了什么内容,来判断是否可以作为“国内外研究现状”来进行详细分析。

27840
来自专栏数据库

《数据库系统概念》15-可扩展动态散列

静态散列要求桶的数目始终固定,那么在确定桶数目和选择散列函数时,如果桶数目过小,随着数据量增加,性能会降低;如果留一定余量,又会带来空间的浪费;或者定期重组散列...

37470
来自专栏简书专栏

Python数据分析及可视化-小测验

本文中测验需要的文件夹下载链接: https://pan.baidu.com/s/1OqFM2TNY75iOST6fBlm6jw 密码: rmbt 下载压缩包...

34220
来自专栏HansBug's Lab

算法模板——Dinic最小费用最大流

实现功能:输入M,N,S,T;接下来M行输入M条弧的信息(包括起点,终点,流量,单位费用);实现功能是求出以S为源点,T为汇点的网络最大流的最小费用 其实相当的...

60960
来自专栏编程札记

数据库内部排序算法之两阶段多路归并排序算法实现

94930

扫码关注云+社区

领取腾讯云代金券