Matlab信號處理基礎

一. 簡介

　　離散傅立葉、離散餘弦和離散小波變換是圖像、音頻信號經常使用基礎操做，時域信號轉換到不一樣變換域之後，會致使不一樣程度的能量集中，信息隱藏利用這個原理在變換域選擇適當位置係數進行修改，嵌入信息，並確保圖像、音頻信號經處理後感官質量無明顯變化。

二. 數學公式

一維離散傅立葉變換對定義

一維離散傅里葉變換：

一維離散傅里葉逆變換：

一維離散餘弦變換對定義

一維離散餘弦正變換：

一維離散餘弦反變換：

 

一維連續小波變換對定義

一維連續小波變換，其總h(t)是小波母函數：

 

一維連續小波逆變換：

二維離散傅立葉變換對定義

二維離散傅里葉變換：

二維離散傅里葉逆變換：

二維離散餘弦變換對定義

二維離散餘弦正變換：

二維離散餘弦反變換：

 

三. 代碼實現

1. 用離散傅里葉變換分析合成音頻和圖像

（1）分析合成音頻文件

第一步：讀取音頻文件數據。

 　　uigetfile() 是文件對話框函數，提供圖形界面供用戶選擇所需文件，返回目標的目錄名和文件名。

 　　函數原型：y= wavread (FILE)

     　　功能：讀取微軟音頻格式（wav）文件內容

    　　輸入參數：file 表示音頻文件名，字符串

    　　返回參數：y 表示音頻樣點，浮點型

第二步：一維離散傅立葉變換

　　fft 函數對輸入參數進行一維離散傅立葉變換並返回其係數，對應頻率從 0 到 fs（採樣頻率），使用 fftshift 將零頻對應係數移至中央。爲了更好地觀察頻譜，可計算離散樣點對應的頻率值。

第三步：一維離散傅立葉逆變換

　　ifft 函數對輸入參數進行一維離散傅立葉逆變換並返回其係數。

第四步：觀察結果

　　figure(n)表示建立第 n 個圖形窗。

　　subplot 是子繪圖函數，第1、二個參數指明子圖像佈局方式，例如，若參數爲 2，3 則表示畫面共分爲 2 行，每行有 3 個子圖像。第三個參數代表子圖像序號，排序順序爲從左至右，從上至下。

　　plot 是繪圖函數，默認使用方式爲 plot(y)，參數 y 是要繪製的數據；若是須要指明圖像橫軸顯示序列，則命令行爲 plot(x, y)，默認方式等同於 plot ([0..len-1], y)，len爲序列y的長度。

例

 

len=40000;
[fn,pn]=uigetfile('*.wav','請選擇音頻文件');
[x,fs]=wavread(strcat(pn,fn),len); %strcat()鏈接生成文件絕對路徑

%一維離散傅里葉變換
xf=fft(x); 
f1=[0:len-1]*fs/len;
xff=fftshift(xf);　　% fftshift()將零頻對應係數移至中央
h1=floor(len/2);　　 %取不大於len/2的最大整數
f2=[-h1:h1]*fs/len; 

%一維離散傅里葉逆變換
xsync=ifft(xf); 
figure;
subplot(2,2,1);plot(x);title('原始圖像');
subplot(2,2,2);plot(xsync);title('synthesiaze audio');
subplot(2,2,3);plot(f1,abs(xf));title('fft coef. of audio'); %abs()求整數的絕對值
subplot(2,2,4);plot(f2(1:len),abs(xff));title('fftshift coef. of audio');

 

 

 

　

（2）分析合成圖像文件

第一步：讀取圖像文件數據

 　　函數原型：A = imread(filename,fmt)

 　　功能：讀取 fmt 指定格式的圖像文件內容

 　　輸入參數：filename 表示圖像文件名，字符串。Fmt 表示圖像文件格式名，字符串，函數支持的圖像格式包括：JPEG，TIFF，GIF，BMP 等等，當參數中不包括文件格式名時，函數嘗試推斷出文件格式。

 　　返回參數：A 表示圖像數據內容，整型。

 　　rgb2gray 函數將 RGB 圖像轉換爲灰度圖。

第二步：二維離散傅立葉變換

　　fft2 函數對輸入參數進行二維離散傅立葉變換並返回其係數，使用 fftshift 將零頻對應係數移至中央。

第三步：二維離散傅立葉逆變換

　　ifft2 函數對輸入參數進行二維離散傅立葉逆變換並返回其係數。

第四步：觀察結果

　　imshow 是二維數據繪圖函數，mesh 經過三維平面顯示數據。

例：

 

[fn,pn]=uigetfile('*.bmp','請選擇圖像文件');
[x,map]=imread(strcat(pn,fn),'bmp');
I=rgb2gray(x);　　%rgb2gray函數將 RGB 圖像轉換爲灰度圖

%二維離散傅立葉變換
xf=fft2(I);
xff=fftshift(xf);

%二維離散傅立葉逆變換
xsync=ifft2(xf);

%結果
figure;
subplot(2,2,1);imshow(x);title('original image');
subplot(2,2,2);imshow(uint8(abs(xsync)));title('synthesize image');
subplot(2,2,3);mesh(abs(xf));title('fft coef. of image');
subplot(2,2,4);mesh(abs(xff));title('fftshift coef. of image');

 

 

 

 

 

2. 用離散餘弦變換分析合成音頻和圖像

（1）分析合成音頻文件

第一步：讀取音頻文件數據。

第二步：一維離散餘弦變換

　　dct 函數對輸入參數進行一維離散餘弦變換並返回其係數，對應頻率從0到 fs（採樣頻率）。

第三步：一維離散餘弦逆變換

　　idct 函數對輸入參數進行一維離散餘弦逆變換並返回其係數。離散餘弦變換經常使用於圖像壓縮，能夠嘗試只使用部分系數重構語言，經過觀察可發現，原始音頻和合成後音頻二者差異不大。

第四步：觀察結果

例：

len=40000;
[fn,pn]=uigetfile('*.wav','請選擇音頻文件');
[x,fs]=wavread(strcat(pn,fn),len);

%一維離散餘弦變換
xf=dct(x);
f1=[0:len-1]*fs/len;
h1=floor(len/2);
f2=[-h1:h1]*fs/len; 

%一維離散餘弦逆變換
xsync=idct(xf);
[row,col]=size(x);
xff=zeros(row,col);
xff(1:row,1:col)=xf(1:row,1:col);
y=idct(xff);

%結果
figure;
subplot(2,2,1);plot(x);title('原始圖像');
subplot(2,2,2);plot(xsync);title('synthesiaze audio');
subplot(2,2,3);plot(f1,abs(xf));title('fft coef. of audio');
subplot(2,2,4);plot(f2(1:len),abs(xff));title('fftshift coef. of audio');

 

 

（2）分析合成圖像文件

第一步：讀取圖像文件數據

第二步：二維離散餘弦變換

　　dct 函數對輸入參數進行二維離散餘弦變換並返回其係數。

第三步：二維離散餘弦逆變換

　　idct2 函數對輸入參數進行二維離散餘弦逆變換並返回其係數。能夠嘗試使用部分系數重構圖像，例如使用係數矩陣中4/5的數據，其它部分置零。

　　爲了保證圖像能正確顯示，使用uint8 對重構圖像原始數據進行了數據類型轉換，確保其取值範圍在 0 到 255 之間。

第四步：觀察結果

    請輸入命令顯示四個子圖，分別是原始圖像、使用所有係數恢復的圖像，使用部分系數恢復的圖像和用三維立體圖方式顯示係數。

例：

[fn,pn]=uigetfile('*.bmp','請選擇圖像文件');
[x,map]=imread(strcat(pn,fn),'bmp');
I=rgb2gray(x);

%二維離散餘弦變換
xf=dct2(I);

%二維離散餘弦逆變換
xsync=uint8(idct2(xf));
[row,col]=size(I);
lenr=round(row*4/5);
lenc=round(col*4/5);
xff=zeros(row,col);
xff(1:lenr,1:lenc)=xf(1:lenr,1:lenc);
y=uint8(idct2(xff));

%結果
figure;
subplot(2,2,1);imshow(x);title('original image');
subplot(2,2,2);imshow(uint8(abs(xsync)));title('synthesize image');
subplot(2,2,3);imshow(uint8(abs(y)));title('part synthesize image');
subplot(2,2,4);mesh(abs(xff));title('fftshift coef. of image');

 

 

 

3. 用離散小波變換分析合成音頻和圖像

（1）分析合成音頻文件

第一步：讀取音頻文件數據。

第二步：一維離散小波變換

　　wavedec 函數對輸入參數進行一維離散小波變換並返回其係數C 和各級係數長度L。第二個參數指明小波變換的級數，第三個參數指明小波變換使用的小波基名稱。

第三步：一維離散小波逆變換

　　waverec 函數對輸入參數進行一維離散小波逆變換並返回其係數。

　　appcoef 函數返回小波係數近似份量，第一個參數 C、第二個參數 L 是 wavedec 的返回參數，爲各級小波係數和其長度，第三個參數指明小波基名稱，第四個參數指明級數。

　　detcoef 函數返回小波係數細節份量，第一個參數 C、第二個參數 L 是 wavedec 的返回參數，爲各級小波係數和其長度，第三個參數指明級數。

第四步：觀察結果

 例：

len=40000;
[fn,pn]=uigetfile('*.wav','請選擇音頻文件');
[x,fs]=wavread(strcat(pn,fn),len);

%一維離散小波變換
[C,L]=wavedec(x,2,'db4');

%一維離散小波逆變換
xsync=waverec(C,L,'db4');
cA2=appcoef(C,L,'db4',2);
cD2=detcoef(C,L,2);
cD1=detcoef(C,L,1);

%結果
figure;
subplot(2,3,1);plot(x);title('原始圖像');
subplot(2,3,2);plot(xsync);title('synthesiaze audio');
subplot(2,3,4);plot(cA2);title('app coef. of audio');
subplot(2,3,5);plot(cD2);title('det coef. of audio');
subplot(2,3,6);plot(cD1);title('det coef. of audio');

 

 

（2）分析合成圖像文件

第一步：讀取圖像文件數據

第二步：二維離散小波變換

　　dwt2函數對輸入參數進行二維一級離散小波變換並返回近似份量，水平細節份量，垂直細節份量和對角線細節份量。若是要對圖像進行多級小波分解，使用wavedec2函數。

第三步：二維離散小波逆變換

　　idwt2 函數對輸入參數進行二維離散小波逆變換並返回其係數。能夠嘗試僅使用近似份量、水平細節份量、垂直細節份量或對角線細節份量重構圖像。

第四步：觀察結果

　　輸入命令顯示六個子圖，分別是原始圖像、使用所有係數恢復的圖像、小波係數近似份量、水平細節份量、垂直細節份量和對角線細節份量。

例：

[fn,pn]=uigetfile('*.bmp','請選擇bmp格式圖像文件');
[x,map]=imread(strcat(pn,fn),'bmp');
I=rgb2gray(x);

%二維離散小波變換
sx=size(I);
[cA1,cH1,cV1,cD1]=dwt2(I,'bior3.7');

%二維離散小波逆變換
xsync=uint8(idwt2(cA1,cH1,cV1,cD1,'bior3.7',sx));
A1=uint8(idwt2(cA1,[],[],[],'bior3.7',sx));
H1=uint8(idwt2([],cH1,[],[],'bior3.7',sx));
V1=uint8(idwt2([],[],cV1,[],'bior3.7',sx));
D1=uint8(idwt2([],[],[],cD1,'bior3.7',sx));

%結果
figure;
subplot(2,3,1);imshow(x);title('original image');
subplot(2,3,2);imshow(uint8(abs(xsync)));title('synthesize image');
subplot(2,3,3);mesh(A1);title('app coef. of image');
subplot(2,3,4);mesh(H1);title('hor coef. of image');
subplot(2,3,5);mesh(V1);title('ver coef. of image');
subplot(2,3,6);mesh(D1);title('dia coef. of image');