實(shí)驗(yàn)三、DFT和DCT及頻域?yàn)V波

1、一 實(shí)驗(yàn)名稱：數(shù)字信號(hào)的 DFT/DCT 及頻域?yàn)V波二 實(shí)驗(yàn)?zāi)康? 熟練掌握數(shù)字信號(hào)（1D）及數(shù)字圖像(2D)離散傅立葉變換（DFT）及離散余弦變換（DCT）方法、基本原理及實(shí)現(xiàn)流程。熟悉兩種變換的性質(zhì)，并能對(duì) DFT 及 DCT 的結(jié)果進(jìn)行必要解釋。2 深入理解離散信號(hào)采樣頻率、奈奎斯特頻率及頻率分辨率等基本概念，弄清它們之間的相互關(guān)系。了解離散傅里葉變換（DFT）中頻率泄露的原因，以及如何盡量減少頻率泄露影響的途徑。3 熟悉和掌握利用 MATLAB 工具進(jìn)行 1D/2D FFT 及 DCT 的基本步驟、MATLAB 函數(shù)使用及對(duì)具體變換的處理流程。4 能熟練應(yīng)用 MATLAB 工具對(duì)數(shù)字圖

2、像進(jìn)行 FFT 及 DCT 處理，并能根據(jù)需要進(jìn)行必要的頻譜分析和可視化顯示。三 實(shí)驗(yàn)原理1、 傅立葉變換l 傅立葉變換：非周期函數(shù)表示為正弦和/或余弦乘以加權(quán)函數(shù)的積分。l 一維連續(xù)Fourier變換對(duì)函數(shù)f(x)進(jìn)行傅立葉變換得到F(u) 逆變換，即將F(u)變換到f(x)為 l 一維離散Fourier變換正變換(DFT) 逆變換(IDFT) l 用幅值和相位表示傅立葉變換 2、 離散余弦變換l 1D-DCT l IDCT變換 l 矩陣形式 四 實(shí)驗(yàn)步驟1. 1D數(shù)字信號(hào)的FFT及頻譜分析給定如下式（1）所示的1D連續(xù)信號(hào)： 1) 設(shè)采樣頻率 =1000Hz，對(duì)信號(hào)進(jìn)行離散化，并畫(huà)出一個(gè)周

3、內(nèi)的信號(hào)振幅隨時(shí)間變化的波形圖。2) 對(duì)離散信號(hào)進(jìn)行傅立葉變換，分別畫(huà)出頻譜中心化及有效頻率范圍（不含負(fù)頻）2種方式下的幅值()隨頻率 變化的分布圖，要求縱橫坐標(biāo)正確標(biāo)注物理量和單位。3) 對(duì)式（1）信號(hào)，加隨機(jī)噪聲，重復(fù)步驟（1）和（2）的處理過(guò)程。4) 通過(guò)對(duì)變換結(jié)果的分析，說(shuō)明采樣頻率、奈奎斯特（Nyquist）頻率（）及采樣時(shí)間間隔T 三者之間的相互關(guān)系，并簡(jiǎn)要描述模擬信號(hào)的采樣定理。圖 1注：加隨機(jī)噪聲的處理只需改變初始信號(hào)即可，其它步驟同理。2. 數(shù)字音頻信號(hào)的DFT1) 讀取一段 0.5s 的預(yù)先錄制的數(shù)字音頻信號(hào)（ “yes.wav”或“no.wav”文件中任選其一），畫(huà)出隨時(shí)

4、間變化的聲波波形圖。2) 對(duì)數(shù)字音頻信號(hào)進(jìn)行離散傅立葉變換（DFT） ，分別畫(huà)出頻譜中心化及有效頻率范圍(不含負(fù)頻)2 種方式下的幅值（）隨頻率變化的分布圖，要求縱橫坐標(biāo)正確標(biāo)注物理量和單位。圖 23. 數(shù)字音頻信號(hào)的 DCT 和 IDCT1) 對(duì)上述音頻信號(hào)做離散余弦變換（DCT），畫(huà)出 DCT 變換系數(shù)（變換結(jié)果）圖，并對(duì)變換結(jié)果進(jìn)行必要的解釋，說(shuō)明 DCT 變換的主要用途。要求按 DCT 原理自行編寫(xiě)實(shí)現(xiàn)代碼，不允許直接調(diào)用 MATLAB 的 dct()函數(shù)。2) 按原理自行一段 MATLAB 代碼，對(duì)第（1）步處理結(jié)果進(jìn)行離散余弦反變換(IDCT)，將計(jì)算結(jié)果與原始音頻信號(hào)進(jìn)行比較，檢

5、驗(yàn)編寫(xiě)代碼的正確性。3) 編寫(xiě)一段 MATLAB 代碼，利用快速傅立葉變換（FFT）程序?qū)崿F(xiàn)快速 DCT 算法（FCT），并將計(jì)算結(jié)果與直接調(diào)用 dct（）的處理結(jié)果進(jìn)行比較，檢驗(yàn)編寫(xiě)代碼的正確性。 圖 34. 綜合應(yīng)用題：實(shí)際信號(hào)的頻譜分析及頻域?yàn)V波1) 編寫(xiě)一從保存在本地磁盤(pán)的文本文件中讀入一實(shí)際數(shù)字信號(hào)，已知該信號(hào)的時(shí)間采樣率為 dt = 2ms。文件中的信號(hào)由 301個(gè)等長(zhǎng)的按列排列的一維列信號(hào)組成，每個(gè)一維列信號(hào)有 251 個(gè)采樣點(diǎn)，信號(hào)實(shí)際計(jì)時(shí)起點(diǎn)為 1800ms，延時(shí)長(zhǎng)度為 L= (251-1)* 2ms =500ms。請(qǐng)讀出其中的某一列信號(hào)，并畫(huà)出該信號(hào)振幅隨時(shí)間變化的波形圖，

6、以 ms 為時(shí)間單位。2) 對(duì)第一步中抽取的其中一列信號(hào)做快速傅里葉變換（FFT），分別畫(huà)出頻譜中心化的對(duì)稱頻譜和只含有正半抽的信號(hào)頻譜圖，并對(duì)該信號(hào)做簡(jiǎn)要的頻譜分析。要求規(guī)范的標(biāo)注縱橫坐標(biāo)實(shí)際物理量和對(duì)應(yīng)的單位。3) 設(shè)定截止頻率 D0=100，試在同一張圖上以不同線型畫(huà)出 n =1，2，4 階下的巴特沃思（ Butterworth）低通濾波器（一維）的頻率響應(yīng)曲線。要求標(biāo)注規(guī)范地縱橫坐標(biāo)實(shí)際物理量和對(duì)應(yīng)的單位。4) 選擇合適的 D0，利用上述 2 階 Butterworth 低通濾波器，對(duì)第（1）步讀取的列信號(hào)進(jìn)行濾波實(shí)驗(yàn)。并分析截止頻率對(duì)濾波效果的影響。圖 4五 實(shí)驗(yàn)結(jié)果及分析圖 51.

7、 1D數(shù)字信號(hào)的FFT及頻譜分析圖 6l 分析：通過(guò)對(duì)變換結(jié)果的分析采樣頻率、奈奎斯特頻率及采樣時(shí)間間隔三者之間有如下的相互關(guān)系 l 模擬信號(hào)的采樣定理：當(dāng)采樣頻率大于連續(xù)信號(hào)最高頻率時(shí)，采樣之后的數(shù)字信號(hào)完整地保留了原始信號(hào)中的信息，是連續(xù)信號(hào)離散化的基本依據(jù)。圖 72. 數(shù)字音頻信號(hào)的DFTl 分析：從上面的頻譜可以看出該音頻的頻率成分主要集中在700Hz以內(nèi)，共有4支峰。圖 83. 數(shù)字音頻信號(hào)的 DCT 和 IDCTl 分析：DCT變換將信號(hào)從時(shí)域轉(zhuǎn)換到變換域上，通過(guò)對(duì)變換后的系數(shù)分析，原能量集中在少數(shù)系數(shù)上，可以提高編碼效率，壓縮數(shù)據(jù)。圖 94. 綜合應(yīng)用題：實(shí)際信號(hào)的頻譜分析及頻域

8、濾波圖 10圖 11圖 12l 該信號(hào)的第二列數(shù)據(jù)的有效頻率主要集中在0-100Hz以內(nèi)l 通過(guò)設(shè)置不同截止頻率的Butterworth低通濾波器，可以看出對(duì)于低頻信號(hào)，截止頻率從100Hz減小，原信號(hào)的有效頻率則被濾去的變多，若是圖像信號(hào)，則表現(xiàn)出圖像變的模糊。六.實(shí)驗(yàn)心得體會(huì)和建議l 心得體會(huì)：通過(guò)這次實(shí)驗(yàn)使我深刻了解了奈奎斯特定理、DFT和DCT的基本原理以及巴特沃斯低通濾波器的構(gòu)造，與此同時(shí)在上機(jī)實(shí)驗(yàn)中熟悉了MATLAB編寫(xiě)FFT和DCT 的基本方法及步驟。在對(duì)相關(guān)頻譜的分析過(guò)程中更加深刻的理解了原理，及相關(guān)用途。l 建議：可以讓大家不調(diào)用FFT函數(shù)，直接編寫(xiě)DCT。七.程序源代碼1.

9、 1D數(shù)字信號(hào)的FFT及頻譜分析fs=1000;N=256;n=0:N-1;t=n/fs;x=2*sin(30*pi*t)+0.5*cos(120*pi*t)+4*sin(240*pi*t);%輸入信號(hào)plot(t,x),xlim(0,1/15);%畫(huà)出一個(gè)周期內(nèi)的信號(hào)振幅隨時(shí)間變化的波形圖%=y=fft(x);mag=abs(y);y=fftshift(y);%頻譜中心化mag0=abs(y);%=M = length(y); % fft頻率軸點(diǎn)數(shù)(may be different)f = (0:M-1)*fs/M; % 頻率采樣序列（矢量）%=fchar = num2str(fs); %

10、采樣率轉(zhuǎn)化為char(本文)nchar = num2str(N); % 樣點(diǎn)數(shù)轉(zhuǎn)化為char(本文)ltext = strcat(fs=,fchar,Hz,N=, nchar, points); %拼title字符%=subplot(311),plot(t,x); % 隨時(shí)間變化的振幅xlabel(t/s);ylabel(振幅);xlim(0,0.1);title(x=2*sin(30*pi*t)+0.5*cos(120*pi*t)+4*sin(240*pi*t);grid on;f0 = f-f(M/2);subplot(312),plot(f0,(mag0)*2/N); % 隨頻率變化的振

11、幅xlabel(頻率/Hz);ylabel(振幅);title(全部頻率：,ltext);grid on;subplot(313),plot(f(1:M/2),(mag(1:M/2)*2/N); % 繪制有效頻譜xlabel(頻率/Hz);ylabel(振幅);title(有效頻率：,ltext);grid on;注：對(duì)于加性噪聲只需將x變?yōu)閤=2*sin(30*pi*t)+0.5*cos(120*pi*t)+4*sin(240*pi*t)+randn(size(t);即可，其它程序同理。 2. 數(shù)字音頻信號(hào)的DFTx,fs=wavread(yes.wav);N=4000;n=0:N-1;t=

12、n/fs;plot(t,x),%畫(huà)出一個(gè)周期內(nèi)的信號(hào)振幅隨時(shí)間變化的波形圖%=y=fft(x);mag=abs(y);y=fftshift(y);%頻譜中心化mag0=abs(y);%=M = length(y); % fft頻率軸點(diǎn)數(shù)(may be different)f = (0:M-1)*fs/M; % 頻率采樣序列（矢量）%=fchar = num2str(fs); % 采樣率轉(zhuǎn)化為char(本文)nchar = num2str(N); % 樣點(diǎn)數(shù)轉(zhuǎn)化為char(本文)ltext = strcat(fs=,fchar,Hz,N=, nchar, points);%拼title字符%=s

13、ubplot(311),plot(t,x); % 隨時(shí)間變化的振幅xlabel(t/s);ylabel(振幅);xlim(0,0.5);title(yes);grid on;f0 = f-f(M/2);subplot(312),plot(f0,(mag0)*2/N); % 隨頻率變化的振幅xlabel(頻率/Hz);ylabel(振幅);title(全部頻率：,ltext);grid on;subplot(313),plot(f(1:M/2),(mag(1:M/2)*2/N); % 繪制有效頻譜xlabel(頻率/Hz);ylabel(振幅);title(有效頻率：,ltext);grid o

14、n;3. 數(shù)字音頻信號(hào)的 DCT 和 IDCTy,Fs= wavread(yes.wav);N=length(y);t=(0:N-1)/Fs;subplot(311),plot(t,y);xlabel(t/s);ylabel(振幅);title(yes波形);G=zeros(N,N);for x=1:N; for u=2:N; G(x,u)=cos(2*(x-1)+1)*(u-1)*pi/(2*N); endendG(1,:)=sqrt(2/N); G(1,1)=sqrt(1/N);F=G*y;subplot(312),plot(1:N)/Fs,F)title(DCT coeffcients)

15、;xlabel(t/s);ylabel(振幅);f=G*F;subplot(313),plot(t,f/Fs*4);title(The recovery of the signal from DCT coefficients)xlabel(t/s);ylabel(振幅);4. 綜合應(yīng)用題：實(shí)際信號(hào)的頻譜分析及頻域?yàn)V波clc,clear,close all;matrix=importdata(seismic_nsamp251_tr301_2ms.txt);p,q=size(matrix);fs=500;r=2;%這里取第2列dt=1/fs;L=(p-1)*dt;%每列的延時(shí)長(zhǎng)度t0=1.8;to

16、_r=t0+L*(r-1);%第r列的開(kāi)始時(shí)間sl=matrix(:,r);%取出第r列的數(shù)據(jù)n=0:p-1;t=to_r+n/fs;% t=n/fs;%=%劃出振幅隨時(shí)間的變化曲線%=figure(1);plot(t,sl); % 隨時(shí)間變化的振幅xlabel(t/ms);ylabel(振幅);title(strcat(第,num2str(r),列);grid on;%=y=fft(sl,p-1);mag=abs(y);y=fftshift(y);%頻譜中心化mag0=abs(y);%=M= length(y); % fft頻率軸點(diǎn)數(shù)(may be different)f = (0:M-1)

17、*fs/M; % 頻率采樣序列（矢量）%=%畫(huà)出中心頻率圖f0 = f-f(M/2);figure;subplot(211),plot(f0,mag0); % 隨頻率變化的振幅xlabel(頻率/Hz);ylabel(振幅);title(全部頻率);grid on;%=subplot(212),plot(f(1:M/2),mag(1:M/2); % 繪制有效頻譜xlabel(頻率/Hz);ylabel(振幅);title(有效頻率);grid on;%=%生成Butterworth矩陣H=zeros(3,M);JS=1 2 4;D0=100;for n=1:M H(1,n)=1/(1+(n-(

18、M)/2)/D0)(2*JS(1); %一階Butterworth矩陣 H(2,n)=1/(1+(n-(M)/2)/D0)(2*JS(2); %二階Butterworth矩陣 H(3,n)=1/(1+(n-(M)/2)/D0)(2*JS(3); %四階Butterworth矩陣end%=%在同一圖中畫(huà)出butterworth函數(shù)figure(3);subplot(311),plot(f0,abs(H(1,:),-b);hold on;subplot(312),plot(f0,abs(H(2,:),-b);hold on;subplot(313),plot(f0,abs(H(3,:),-b);hold on;xlim(f0(1)