最近著手把CSK移植到DSP中，先看一些DSP中圖像處理的一些例子，第一件事當然就是怎么把圖像數(shù)據(jù)倒入CCS工程中了，去年倒是用過一點CCS，再拿起來已經(jīng)忘得差不多了，這篇文章主要記錄一些學習的過程：
目錄
- DSP導入圖像數(shù)據(jù)
- 窗函數(shù)和加窗
- 定點數(shù)和浮點數(shù)
- 圖像的FFT運算
- 圖像的IFFT運算
一. DSP導入圖像數(shù)據(jù)
搞了一下午大概可以了，主要是如何導入數(shù)據(jù)，如何利用CCS的Image analyzer來做顯示。
1. 準備dat形式的數(shù)據(jù)
CCS導入數(shù)據(jù)是dat格式的，利用matlab把數(shù)據(jù)保存到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的格式規(guī)則，要用到的話再細究。
2. 準備內(nèi)存，導入數(shù)據(jù)
可以利用CCS的load memory功能，這個功能是在調(diào)試的時候可用的，所以先開辟內(nèi)存，然后再load memory。這里也比較簡單：
點擊next選擇內(nèi)存位置和長度：
注意這里數(shù)組不能太長，可能會造成溢出，我把堆棧設置改成0xFFFF，這里用2500的長度就可以了，具體這塊還沒有搞清楚，50*50的圖片也夠我測試了。點擊finish就可以完成導入了，然后可以在memory browser里查看數(shù)據(jù)：
查看圖像。
查看圖像用tools->image analyzer這個工具，會彈出兩個框，一個屬性框properties設置一些參數(shù)：
然后在對應的圖像窗口點refresh就可以看到圖像了：
加載圖像加載進去倒是可以，但是是訪問不了的，寫了個很簡單的閾值處理跑不了。
替代方案
最后也沒有發(fā)現(xiàn)到底是什么樣的原因，但是在李老師（李翔宇）的指導下，把圖片數(shù)據(jù)的首地址給定，最后給在的是0x00850000，這樣就可以了，老師說也 不知道什么原因，先放下吧，最起碼這樣可以進行一些簡單的算法仿真了。
具體這樣：
unsigned char *img=(unsigned char *)0x0850000;
再load memory就沒有問題了，也可以進行操作。
二.窗函數(shù)的實現(xiàn)和導入
CSK的實現(xiàn)過程中要用到兩種窗函數(shù)，分別是高斯和漢明，這兩種窗函數(shù)可以利用matlab提前生成好，然后作為頭文件來導入到CCS工程中。這個實現(xiàn)起來也不難。
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建立的頭文件中了，數(shù)據(jù)類型選擇float類型。我試著導入了一個漢明窗，64*64的，拉成一維數(shù)組后導入是這個樣子：
CCS中如果想看二維圖像長什么樣子的話還得把float轉(zhuǎn)換成8位int，挺麻煩的，所以就這樣看看吧。
然后我們把窗函數(shù)加在圖像上看一看長什么樣子，我這里換了一副圖像，所有的操作都是針對6464的圖像來的。
導入的圖像為
加窗是要把兩個數(shù)組相乘，直接用循環(huán)把兩個數(shù)組乘起來就可以了（element-wise）,floatuchar=float，這樣得到的結果如下：
然后就是把乘的結果轉(zhuǎn)換為uchar型的來顯示，首先定義一個存放轉(zhuǎn)換之后結果的數(shù)組，然后用循環(huán)逐一轉(zhuǎn)化（這種應該都是可以用多核進行優(yōu)化并行計算的）我一開始是這么寫的：
for(i=0;i<4096;i++)
{
*(img_smooth++)=*(img_with_window++); //轉(zhuǎn)換為uchar
}
會莫名奇妙出現(xiàn)問題，不僅img_smooth的結果不對，還會把原先img_with_window的數(shù)據(jù)給沖掉。數(shù)組地址我都是手動給定的，保證不會沖突的，也沒找到什么原因。（如有人知道還望指點！）
img_smooth就變成這樣了，這肯定是不對的。而且img_with_window就全變成0了。
后來換了一種方法用下標來做，就沒有問題了，真是神奇。
轉(zhuǎn)換為uchar之后長這樣，對比可以知道這個是正確的，這個時候就可以看加窗之后的圖像了：
都是沒有問題的，這個smooth的命名的原因就是因為上次寫均值濾波時用的那個變量，還沒有改。
三.定點數(shù)和浮點數(shù)的區(qū)別
PC編程很少遇到這么細節(jié)的問題，但是DSP上就不同了，以前只知道定點數(shù)需要定標，浮點數(shù)是采用類似于科學計數(shù)法的一種方法，具體的細節(jié)就不清楚了，DSP還有定點和浮點之分，所以把這里的細節(jié)看了看，然后整理下。
定點數(shù)
所謂定點數(shù)就是約定小數(shù)點的位置是不變的，這樣通常將定點數(shù)據(jù)表示成純小數(shù)或者純整數(shù)，純小數(shù)的話，就把小數(shù)點固定在數(shù)值部分的最高位之前；純小數(shù)的話，則把小數(shù)點固定在數(shù)值部分的最后面，如圖：
小數(shù)點在機器中是不表示出來的，是一個提前約定好的位置，對于一臺計算機，一旦確定了小數(shù)點的位置，就不再改變。
必須明確的一點就是計算機是無法計算小數(shù)的，所有的小數(shù)計算都是通整數(shù)計算完成的，這就導致事先約定小數(shù)點的位置尤為重要，這就是定標，在定點運算中，定標很重要。
為了簡單，我們用16位的來說吧，也就是說參與運算的數(shù)都是16位的整型數(shù)，那么如何處理小數(shù)的關鍵就在于如何確定一個數(shù)據(jù)小數(shù)點的位置，這樣的話數(shù)據(jù)的精度和數(shù)據(jù)的范圍就是一對矛盾了，精度越高，范圍則會越小。
數(shù)的標定有Q和S兩種表示方法：
對于16位的來說，最高位為符號保留位，那么根據(jù)定標的不同，有Q0到Q15共16中定標方式，精度依次增高，范圍依次降低。
例：
0010 0000 0000 0000b 表示8192 Q0定標
而同樣的一個數(shù)：Q15定標
0010 0000 0000 0000b 表示0.25
這個圖片不是特別清楚，可以看下：
DSPlib中的fft32x32函數(shù)就是以Q15定標的，所以如果要用這個函數(shù)，就要對數(shù)據(jù)進行轉(zhuǎn)換。這個可以自己手動定標，也可以用DSPlib里的函數(shù)，我找到一個函數(shù)叫做DSP_fltoq15（）,這個是float到Q15的一個轉(zhuǎn)換函數(shù)，是可以進行這樣的轉(zhuǎn)換的。這個在CCSdsplib的文檔中就有。
這個算法是直接給出的，也很簡單，實際上是一個移位（乘法就相當于移位），注意調(diào)用之前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;
}
}
至于浮點數(shù)怎么判斷溢出和運算規(guī)則就暫且不細究了，我這次的應用只涉及到所用的算法是要求Q15格式，所以需要手動轉(zhuǎn)換。
浮點數(shù)
正是由于定點數(shù)的表示方法太過僵硬，固定的小數(shù)點位置決定了固定位數(shù)的整數(shù)部分和小數(shù)部分，不利于表示特別大和特別小的數(shù)，現(xiàn)在絕大多數(shù)的現(xiàn)代計算機系統(tǒng)都采納了浮點數(shù)表達實數(shù)。
浮點數(shù)用一個尾數(shù)，一個基數(shù)，一個指數(shù)來表示一個數(shù)。
比如123.45用十進制的科學技術法可以表示為：1.2345*10^2。其中1.2345就是尾數(shù)，10是基數(shù)，2是指數(shù)，這樣可以更加靈活得表示更大范圍的實數(shù)。
IEEE浮點數(shù)
IEEE標準規(guī)定了兩種基本的浮點格式：單精度和雙精度，其中單精度具有23位有效數(shù)字，總共占用32位，那么自然是有8位指數(shù)，1位是保留符號位，雙精度共64位，其中有52位尾數(shù)，11位指數(shù)，1位是保留符號位。
現(xiàn)代的X86架構早已摒棄定點數(shù)，在可以進行浮點數(shù)進行計算的計算環(huán)境中，可以自由得利用float表示小數(shù)進行計算。
細節(jié)不說了，一般不是很涉及底層的話也用不到這些知識。
四.圖像的FFT運算
我在另一篇文章里總結了，二維的FFT可以轉(zhuǎn)換為一維的FFT進行計算，即先對每一列進行FFT運算，然后再第一次FFT運算的結果上進行列FFT變換，具體的證明去看那篇文章，這樣二維的就可以用一維的來做，利用的是dsplib的fft32x32函數(shù)，這個函數(shù)可以快速計算FFT，先來看這個函數(shù)：
這里的w是提前生成的，利用dsolib里的一個exe文件就可以了。
exe
用cmd來打開這個：
直接看截圖吧，這個是一位老師指導的，使用說明也給的很清楚，這樣之后就會生成一個文件，這個文件可以是頭文件，也可以是c文件，我是直接生成c文件，把系數(shù)復制到我的程序里了。
另外x[2nx]和y[2nx]都是32位的數(shù)據(jù)，虛部緊跟實部連續(xù)存儲，32位數(shù)據(jù)采用Q15定點，所以這里涉及到一個數(shù)據(jù)轉(zhuǎn)換，整個fft二維算法的大框架我是借鑒我們組另外一個老師的，所以不便放出來，需要把數(shù)據(jù)轉(zhuǎn)換成Q15格式的，這里的轉(zhuǎn)換建議先歸一化到[-1.1),然后左移15位賦值給int型的，這樣就能得到Q15的數(shù)據(jù)。
nx是FFT變換的點數(shù)，無需多說。
看下結果：對于一個64*64的圖像，我用matlab和ccs分別進行了計算，先歸一化，然后去做FFT，這里只看實部：
精度是有損失的，但是看著前幾十個數(shù)的結果還算對的，但是我把這個結果導入到matlab里和matlab計算的fft的結果相比（我是把數(shù)據(jù)拉長做了個差），發(fā)現(xiàn)一個問題：
實際上是每四行的第一行結果是正確的，其他行的結果是有問題的，這個我估計是它的頻譜是亂序的，具體什么原因現(xiàn)在還不清楚，等我下午再去研究研究，應該是順序出了問題。
我先把中間的結果輸出，即每一列先做FFT的結果保存起來看一下，我發(fā)現(xiàn)這個結果是在正確的，在matlab里和上面同樣做了差，兩位小數(shù)精度下是完全準確的。
所以說這個fft32x32的函數(shù)是不要求輸入亂序的，這個我翻手冊的時候也沒發(fā)現(xiàn)要求輸入亂序。我試了試暫時還沒有找到結果錯誤的原因。我也試了下看是不是亂序，以matlab_fft變換的第二列為標準，在ccs的結果中找最相似的一整列，實驗證明并不存在這樣的一整列。說明不是亂序的結果。
搞了半個下午，還是沒有發(fā)現(xiàn)問題所在，我去干干別的再來搞這個，是在不行要去請教下老師！
18年1月19號上午更新：
剛才正在看別的東西，突然靈機一動，想了想是不是內(nèi)存地址互相沖突呢？于是趕緊打開CCS試了一下，第一次沒成功，趕緊檢查發(fā)現(xiàn)還是有一個變量的地址沒有完全指定，導致兩個變量的地址是完全一樣的，我的程序里是xtemp和ytemp,指定了以后趕緊去看第二行第一個數(shù)，果然就和以前的不同了，導入matlab簡單畫了個圖，perfect！
matlab和CCS傅里葉變換
誤差基本在兩個三個小數(shù)點之后，可以接受！
五.圖像的IFFT運算
FFT成功之后，還是花了半個下午的時間把測試的代碼整理成函數(shù)放到頭文件里，因為代碼不多，就沒有新建源文件，之所以花了半個下午是因為當時一個內(nèi)存死活寫不進去數(shù)，對于一個int型的數(shù)據(jù)取float直接就是零，搞得我都想把電腦砸了，后來重啟了一下CCS竟然好了，神奇的錯誤，這也是現(xiàn)在不想做太多底層的原因，特別是和硬件相關的。
這兩天爸媽來西安了，中午剛送走，從車站趕回來剛趕上下午上班時間，實際上都困得不行了，就想去睡覺，到了辦公室把早上改的IFFT看了一下，其實IFFT就只需要把FFT的代碼中的dsp_fft32x32改成dsp_ifft32x32就可以了，其他的東西都不用改，因為其他得輔助代碼只是為了先行后列這樣操作罷了。
但是算的結果確差了很多，早上沒想到底是怎么回事，下午回來把數(shù)據(jù)導入到matlab里reshape顯示一看，竟然是原圖，那么肯定是數(shù)據(jù)整個擴大了一個整倍數(shù)，找了一組對應的數(shù)字一除，4095.6，簡直完美，這才想起逆傅里葉變換是有個系數(shù)的，DSPlib中的ifft函數(shù)并沒有包含這個系數(shù)，因為是行列做了兩次，所以應該是1/64*64=1/4096，于是在最后除上這個數(shù)就ok了。
這樣這個CSK算法最難的地方就過去了。
2018/1/24更新
六：C674Xdsplib導入及使用
昨天碰見老師了，問我做的怎么樣了，交流了一下決定先用浮點的DSP來做，如果速度可以的話再移植到定點上去，所以前面用定點做的要改一改，于是又去下載了C674x的daplib的庫，這里面有一個dsp的函數(shù)是浮點的，還有一些矩陣乘法數(shù)組乘法之類的東西，昨天一天都在鼓搗這個東西，因為網(wǎng)上的資料實在是太少，搞得有一陣子都想把電腦砸了，不過趕在晚飯之前還是弄好了，不過就是這樣也很氣，感覺做了許多無用功，沒人帶做東西真的累。一晚上做夢也都是浮點定點這點破事，一天好像都在看內(nèi)存里存的什么東西。
不過抱怨歸抱怨，還是把這里總結下，如果有人能看到且要用到，希望有用。
全部復制過去，然后項目里吧include文件夾import進去，用的時候要包含相應的頭文件，這個現(xiàn)在不用了不說了，很簡單。
重點說下dspC674x的dsplib的調(diào)用，這個方法可以用在許多l(xiāng)ib的導入上，比如mathlib，DSPc674的dsplib需要用到mathlib，所以這兩個都導入了，首先去下載兩個庫，TI的官網(wǎng)上就可以直接下，exe格式的，安裝完成后是一個文件夾，長這樣：
readme里說了這些東西都是干嘛的，主要的東西都在packages是里，這個里面是：
CCS，項目的properity里：
把這幾個都添加進去。
還有一個：
linker里面把lib添加進去，就可以用了，用的時候只需要添加#include"dsplib.h"就可以了，這個頭文件打開看一下就明白了，它實際上是一個鏈接頭文件，包含了dsplib里所有的函數(shù)，當然也可以用哪個包含哪個。
這樣準備工作就完成了，這個浮點庫里做FFT用的是這樣一個函數(shù)：
點開可以看原圖，輸入輸出都是float形式的，虛部緊跟實部保存，N是數(shù)目，這里面主要有兩個參數(shù)需要解釋，第一個是brev是一個亂序表，第二個是ptr_w是旋轉(zhuǎn)矩陣。
brev用這個就可以了，也有說這個參數(shù)根本就沒用，我還是寫上了。
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
};
另外有一個函數(shù)：
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;
}
}
}
我想直接調(diào)用這個函數(shù)出不來，覺得還不如自己生成數(shù)據(jù)拷貝進來方便，F(xiàn)FT和IFFT用的是一套w，所以只需要生成一次復制進自己的頭文件或者源文件就行了。我是直接在VS里把這個函數(shù)復制進去，生成的數(shù)據(jù)拷貝到CCS的頭文件中，竟然出奇得順利就完成了。
未完待續(xù)---
2018.7.16： 這個暫時就終結了，這個算法最終大部分的函數(shù)都寫完了但是沒有調(diào)通，主要的原因可能還是在我比較抗拒DSP吧，所以做的成果都打包發(fā)給老師了，老師讓一個職工往下做了，希望可以有些不錯的結果。
