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

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

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

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

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

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

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

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

9、 1D數字信號的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);%輸入信號plot(t,x),xlim(0,1/15);%畫出一個周期內的信號振幅隨時間變化的波形圖%=y=fft(x);mag=abs(y);y=fftshift(y);%頻譜中心化mag0=abs(y);%=M = length(y); % fft頻率軸點數(may be different)f = (0:M-1)*fs/M; % 頻率采樣序列（矢量）%=fchar = num2str(fs); %

10、采樣率轉化為char(本文)nchar = num2str(N); % 樣點數轉化為char(本文)ltext = strcat(fs=,fchar,Hz,N=, nchar, points); %拼title字符%=subplot(311),plot(t,x); % 隨時間變化的振幅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;注：對于加性噪聲只需將x變為x=2*sin(30*pi*t)+0.5*cos(120*pi*t)+4*sin(240*pi*t)+randn(size(t);即可，其它程序同理。 2. 數字音頻信號的DFTx,fs=wavread(yes.wav);N=4000;n=0:N-1;t=

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

13、ubplot(311),plot(t,x); % 隨時間變化的振幅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. 數字音頻信號的 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. 綜合應用題：實際信號的頻譜分析及頻域濾波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;%每列的延時長度t0=1.8;to

16、_r=t0+L*(r-1);%第r列的開始時間sl=matrix(:,r);%取出第r列的數據n=0:p-1;t=to_r+n/fs;% t=n/fs;%=%劃出振幅隨時間的變化曲線%=figure(1);plot(t,sl); % 隨時間變化的振幅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頻率軸點數(may be different)f = (0:M-1)

17、*fs/M; % 頻率采樣序列（矢量）%=%畫出中心頻率圖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%=%在同一圖中畫出butterworth函數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)