分析與處理實(shí)驗(yàn):Matlab波形到頻域分析的完整實(shí)操記錄)
簡介這份資源是南京郵電大學(xué)「信號(hào)分析與處理實(shí)驗(yàn)」課程的完整實(shí)驗(yàn)報(bào)告面向正在修讀數(shù)字信號(hào)處理、信號(hào)與系統(tǒng)相關(guān)課程的高校學(xué)生以及需要借助 MATLAB 完成信號(hào)處理實(shí)驗(yàn)與課程設(shè)計(jì)的自學(xué)者。報(bào)告圍繞信號(hào)的產(chǎn)生和運(yùn)算、連續(xù)時(shí)間信號(hào)的頻域分析、信號(hào)的時(shí)域采樣和重建、離散傅里葉變換與應(yīng)用等核心實(shí)驗(yàn)展開并延伸至連續(xù)時(shí)間系統(tǒng)的頻域與復(fù)頻域分析、數(shù)字濾波器設(shè)計(jì)及信號(hào)處理綜合應(yīng)用配有可運(yùn)行的 MATLAB 程序示例與圖形結(jié)果便于對(duì)照理解時(shí)域與頻域特性。資源包共 1 個(gè) docx 文件約 1.02MB內(nèi)容涵蓋實(shí)驗(yàn)?zāi)康摹?shí)驗(yàn)內(nèi)容、程序代碼與結(jié)果分析結(jié)構(gòu)完整、條理清晰。目前已有 952 人學(xué)習(xí)下載適合作為實(shí)驗(yàn)預(yù)習(xí)、報(bào)告撰寫參考與 DSP 入門練習(xí)的實(shí)用材料。1. 信號(hào)分析與處理實(shí)驗(yàn)從Matlab波形到頻域分析的完整實(shí)操記錄很多同學(xué)學(xué)完《信號(hào)與系統(tǒng)》理論課面對(duì)傅里葉級(jí)數(shù)、采樣定理、DFT這些概念還是一頭霧水公式推得再熟一到實(shí)際信號(hào)上就不知道從哪下手。這份實(shí)驗(yàn)報(bào)告記錄的就是把抽象公式變成Matlab波形的過程——從最基礎(chǔ)的階躍序列、復(fù)指數(shù)序列生成到周期矩形脈沖的頻譜分析再到時(shí)域采樣與信號(hào)重建最后落到數(shù)字濾波器設(shè)計(jì)。四個(gè)實(shí)驗(yàn)層層遞進(jìn)每個(gè)都配有可直接運(yùn)行的Matlab代碼和結(jié)果分析。如果你正在學(xué)數(shù)字信號(hào)處理DSP這門課或者需要一份能照著跑通的實(shí)驗(yàn)參考這份材料能幫你把理論和代碼之間的那層窗戶紙捅破。2. 實(shí)驗(yàn)環(huán)境與項(xiàng)目結(jié)構(gòu)先搞清楚手里有什么2.1 運(yùn)行環(huán)境與依賴這套實(shí)驗(yàn)全部基于Matlab完成不依賴任何額外工具箱。我實(shí)測(cè)在Matlab R2018b及以上版本都能直接跑低版本理論上也沒問題因?yàn)橛玫降亩际腔A(chǔ)函數(shù)stem、plot、subplot、conv、fft、sinc、filter這些。不需要Simulink不需要DSP System Toolbox純腳本就能搞定。如果你用的是Octave大部分代碼也能跑但freqz和butter這類濾波器設(shè)計(jì)函數(shù)可能需要安裝signal包。我一般建議直接用Matlab省得折騰兼容性問題。2.2 實(shí)驗(yàn)?zāi)K劃分整份報(bào)告按實(shí)驗(yàn)序號(hào)組織每個(gè)實(shí)驗(yàn)對(duì)應(yīng)一個(gè)獨(dú)立文件夾里面包含腳本文件和結(jié)果圖。目錄結(jié)構(gòu)大致如下信號(hào)分析與處理實(shí)驗(yàn)/ ├── 實(shí)驗(yàn)一_信號(hào)的產(chǎn)生和運(yùn)算/ │ ├── P1_1.m % 單位階躍序列 │ ├── P1_2.m % 復(fù)指數(shù)序列 │ ├── P1_3.m % 實(shí)指數(shù)序列 │ ├── P1_4.m % 正弦序列 │ ├── P1_5.m % 序列加乘運(yùn)算 │ ├── P1_6.m % 奇偶分量分解 │ └── P1_7.m % 線性卷積 ├── 實(shí)驗(yàn)二_連續(xù)時(shí)間信號(hào)頻域分析/ │ ├── P2_1.m % 周期矩形脈沖傅里葉級(jí)數(shù) │ ├── P2_2.m % 諧波疊加 │ └── P2_4.m % 吉布斯現(xiàn)象過沖計(jì)算 ├── 實(shí)驗(yàn)三_時(shí)域采樣與重建/ │ ├── P3_1.m % 采樣與頻譜分析 │ └── P3_2.m % 信號(hào)重建 └── 實(shí)驗(yàn)四_離散傅里葉變換與應(yīng)用/ ├── P4_1.m % DFT實(shí)現(xiàn) └── P4_2.m % 濾波器設(shè)計(jì)每個(gè)腳本都是獨(dú)立可運(yùn)行的變量命名清晰沒有跨文件依賴。這種組織方式對(duì)新手很友好——你不需要理解整個(gè)工程打開哪個(gè)文件就跑哪個(gè)。2.3 實(shí)驗(yàn)報(bào)告的核心產(chǎn)出每個(gè)實(shí)驗(yàn)最終要交付三樣?xùn)|西一是可運(yùn)行的Matlab代碼二是波形圖或頻譜圖三是對(duì)結(jié)果的分析說明。報(bào)告里對(duì)每個(gè)問題的分析都寫得比較具體比如采樣定理那部分不僅畫出了混疊前后的頻譜對(duì)比還解釋了為什么Ts2時(shí)會(huì)發(fā)生混疊而Tspi/2時(shí)不會(huì)。這種“代碼圖解釋”的結(jié)構(gòu)拿來當(dāng)實(shí)驗(yàn)?zāi)0寤蛘咦詫W(xué)材料都合適。3. 核心實(shí)驗(yàn)拆解從信號(hào)生成到頻域分析的代碼實(shí)現(xiàn)3.1 離散序列的生成與基本運(yùn)算實(shí)驗(yàn)一的核心是讓你熟悉Matlab里怎么表示離散時(shí)間信號(hào)。離散序列和連續(xù)信號(hào)最大的區(qū)別在于你得同時(shí)定義時(shí)間軸向量和幅度向量。比如要畫單位階躍序列u[n]代碼是這樣的% 產(chǎn)生單位階躍序列 u[n] 和 u[n7] clear; close all; n -10:20; % 時(shí)間軸范圍 x [zeros(1,10) 1 ones(1,20)]; % u[n]: n0時(shí)為0, n0時(shí)為1 subplot(1,2,1); stem(n, x); xlabel(n); ylabel(u[n]); title(單位階躍序列 u[n]); M 7; % 移位量 xd [zeros(1,10M) 1 ones(1,20-M)]; % u[n7]: 左移7位 subplot(1,2,2); stem(n, xd); xlabel(n); ylabel(u[n7]); title(移位后的階躍序列 u[n7]);這里的關(guān)鍵點(diǎn)在于zeros和ones的拼接方式。n-10:20一共31個(gè)點(diǎn)zeros(1,10)對(duì)應(yīng)n-10到n-11對(duì)應(yīng)n0ones(1,20)對(duì)應(yīng)n1到n20。移位操作u[n7]本質(zhì)上是把序列左移7位所以前面補(bǔ)7個(gè)零后面少7個(gè)一。這個(gè)邏輯搞清楚了后面所有序列的移位操作都是一個(gè)套路。復(fù)指數(shù)序列稍微繞一點(diǎn)因?yàn)樯婕皬?fù)數(shù)運(yùn)算% 復(fù)指數(shù)序列 x[n] 1.5 * exp((0.2 j3)n), 0 n 15 clear; close all; c 0.2 3*1i; % 復(fù)數(shù)指數(shù)系數(shù) K 1.5; % 幅度系數(shù) n 0:15; x K * exp(c * n); % 生成復(fù)指數(shù)序列 subplot(1,2,1); stem(n, real(x)); % 取實(shí)部 xlabel(n); ylabel(幅度); title(實(shí)部); subplot(1,2,2); stem(n, imag(x)); % 取虛部 xlabel(n); ylabel(幅度); title(虛部);real(x)和imag(x)分別提取實(shí)部和虛部。復(fù)指數(shù)序列在DSP里非常重要因?yàn)樗荄FT的基函數(shù)。這里c 0.2 3i實(shí)部0.2控制包絡(luò)的衰減或增長虛部3控制振蕩頻率。你可以試著改一下c的值觀察波形怎么變——比如把虛部改成0.5振蕩會(huì)明顯變慢。3.2 序列的卷積與奇偶分解卷積是信號(hào)處理里最核心的運(yùn)算之一。實(shí)驗(yàn)一里用conv函數(shù)計(jì)算了兩個(gè)有限長序列的線性卷積% 計(jì)算 x[n] 和 h[n] 的線性卷積 clear; close all; nx 0:14; x (-1).^nx; % x[n] (-1)^n, 0n14 nh 0:7; w 0.25 * pi; h cos(w * nh); % h[n] cos(0.25*pi*n), 0n7 y conv(x, h); % 線性卷積 n 0:21; % 卷積結(jié)果長度 15 8 - 1 22 subplot(3,1,1); stem(nx, x); title(x[n]); subplot(3,1,2); stem(nh, h); title(h[n]); subplot(3,1,3); stem(n, y); title(y[n] x[n] * h[n]);conv函數(shù)的輸出長度是length(x) length(h) - 1這里是158-122個(gè)點(diǎn)。注意時(shí)間軸n要重新定義不能直接用nx或nh。這個(gè)細(xì)節(jié)很容易翻車——我見過不少同學(xué)畫卷積結(jié)果時(shí)時(shí)間軸對(duì)不上波形整體偏移了。奇偶分解的代碼也值得看一下% 序列的奇偶分量分解 clear; close all; n -2:2; x [1 2 4 -5 -2]; y fliplr(x); % 時(shí)間反轉(zhuǎn) xe 0.5 * (x y); % 偶部 xo 0.5 * (x - y); % 奇部 subplot(3,1,1); stem(n, x); title(x[n]); subplot(3,1,2); stem(n, xe); title(偶部 x_e[n]); subplot(3,1,3); stem(n, xo); title(奇部 x_o[n]);fliplr是左右翻轉(zhuǎn)對(duì)應(yīng)時(shí)間反轉(zhuǎn)x[-n]。偶部是(x[n]x[-n])/2奇部是(x[n]-x[-n])/2。這個(gè)分解在分析系統(tǒng)對(duì)稱性時(shí)很有用比如實(shí)偶信號(hào)的傅里葉變換也是實(shí)偶函數(shù)。3.3 周期矩形脈沖的傅里葉級(jí)數(shù)與吉布斯現(xiàn)象實(shí)驗(yàn)二的重頭戲是周期矩形脈沖的頻譜分析。指數(shù)形式的傅里葉級(jí)數(shù)系數(shù)公式是$$c_k \frac{A\tau}{T} \cdot \text{sinc}\left(\frac{k\Omega_0\tau}{2\pi}\right)$$代碼實(shí)現(xiàn)時(shí)用sinc函數(shù)直接算% 周期矩形脈沖的傅里葉級(jí)數(shù)系數(shù) A 1; tau 0.1; T 0.5; Omega0 2*pi/T; k -20:20; ck A*tau/T * sinc(k*Omega0*tau/(2*pi)); subplot(2,1,1); stem(k, abs(ck)); title(雙邊幅度譜); subplot(2,1,2); stem(k, angle(ck)); title(相位譜);這里k取-20到20覆蓋了主要諧波分量。sinc函數(shù)在Matlab里定義為sin(pi*x)/(pi*x)和信號(hào)處理里的歸一化sinc一致。如果你把T從0.5改成1或2會(huì)發(fā)現(xiàn)譜線間隔變密、幅度降低——這正是報(bào)告里問題2.2要你觀察的現(xiàn)象。吉布斯現(xiàn)象的計(jì)算用了一個(gè)循環(huán)來累加諧波function p2_4(N) t -2:1e-4:2; A 1; tau 1; T 2; Omega0 2*pi/T; c0 A*tau/T; f_N c0 * ones(1, length(t)); for k 1:N f_N f_N 2*A*tau/T * sinc(k*Omega0*tau/(2*pi)) * cos(k*Omega0*t); end overshoot max(f_N); fprintf(The maximum overshoot is %.4f\n, overshoot); end當(dāng)N增大時(shí)過沖峰值趨近于一個(gè)常數(shù)大約在跳變幅度的9%左右。這個(gè)函數(shù)直接返回過沖最大值方便你驗(yàn)證不同N下的收斂情況。3.4 采樣定理驗(yàn)證與信號(hào)重建實(shí)驗(yàn)三用升余弦脈沖信號(hào)來驗(yàn)證采樣定理。核心邏輯是先定義連續(xù)信號(hào)的頻譜然后以不同采樣間隔Ts采樣觀察頻譜是否混疊。% 采樣與頻譜分析 Ts pi/2; % 采樣間隔 t -10:Ts:10; % 升余弦脈沖信號(hào)生成 x ...; % 具體表達(dá)式見報(bào)告 % 計(jì)算并繪制頻譜 X fftshift(fft(x)); f (-length(x)/2:length(x)/2-1) / (length(x)*Ts); plot(f, abs(X));當(dāng)Ts pi/2時(shí)采樣頻率滿足奈奎斯特條件頻譜沒有混疊當(dāng)Ts 2時(shí)采樣頻率低于兩倍截止頻率頻譜出現(xiàn)混疊。報(bào)告里對(duì)這兩個(gè)情況的對(duì)比圖和分析都很清楚。信號(hào)重建部分用低通濾波器從采樣信號(hào)中恢復(fù)原信號(hào)% 信號(hào)重建 Wc 2; % 低通濾波器截止頻率 % 設(shè)計(jì)理想低通濾波器 h Wc/pi * sinc(Wc/pi * (n - delay)); % 卷積實(shí)現(xiàn)重建 xr conv(xs, h);重建誤差的大小直接反映了采樣定理是否被滿足。Ts2時(shí)誤差明顯增大因?yàn)榛殳B導(dǎo)致高頻分量無法恢復(fù)。4. 避坑與常見問題跑代碼時(shí)容易翻車的幾個(gè)地方4.1 時(shí)間軸對(duì)不上導(dǎo)致波形偏移現(xiàn)象畫卷積結(jié)果或移位序列時(shí)波形整體左移或右移和理論位置對(duì)不上。原因stem或plot的橫軸向量和信號(hào)向量長度不一致或者移位操作后沒有重新定義時(shí)間軸。解決每次生成新序列后先確認(rèn)length(n) length(x)。卷積后用n 0:length(y)-1重新定義時(shí)間軸。移位操作時(shí)時(shí)間軸不變只改幅度向量的補(bǔ)零位置。4.2 復(fù)數(shù)序列取實(shí)部虛部時(shí)忘記用real/imag現(xiàn)象stem畫復(fù)指數(shù)序列時(shí)報(bào)錯(cuò)或畫出奇怪的結(jié)果。原因Matlab的stem不能直接畫復(fù)數(shù)必須顯式取實(shí)部或虛部。解決用real(x)和imag(x)分別提取。如果想畫幅度和相位用abs(x)和angle(x)。4.3 采樣頻率設(shè)置不當(dāng)導(dǎo)致頻譜混疊現(xiàn)象頻譜圖上出現(xiàn)明顯的頻率折疊高頻分量跑到低頻區(qū)域。原因采樣間隔Ts大于奈奎斯特間隔不滿足采樣定理。解決先計(jì)算信號(hào)的最高頻率fm然后確保Ts 1/(2*fm)。報(bào)告里升余弦脈沖的截止頻率是Wm2所以奈奎斯特間隔是pi/2。如果你不確定信號(hào)的帶寬先把Ts設(shè)小一點(diǎn)逐步增大觀察頻譜變化。4.4 fft后的頻率軸計(jì)算錯(cuò)誤現(xiàn)象頻譜圖的橫軸頻率值不對(duì)峰值位置和理論值偏差很大。原因fft輸出的頻率分辨率是fs/N其中fs1/TsN是采樣點(diǎn)數(shù)。頻率軸應(yīng)該從-fs/2到fs/2。解決用f (-N/2:N/2-1) * (fs/N)計(jì)算頻率軸配合fftshift使用。如果只畫正頻率部分取前N/21個(gè)點(diǎn)頻率軸從0到fs/2。4.5 濾波器設(shè)計(jì)時(shí)階數(shù)選擇不當(dāng)現(xiàn)象重建信號(hào)要么過渡帶太寬要么通帶紋波太大。原因?yàn)V波器階數(shù)太低導(dǎo)致過渡帶不夠陡峭階數(shù)太高又可能引入數(shù)值不穩(wěn)定。解決先用buttord或cheb1ord估算所需階數(shù)再用butter或cheby1設(shè)計(jì)。對(duì)于理想低通重建直接用sinc函數(shù)生成沖激響應(yīng)但要注意截?cái)嚅L度——截?cái)嗵虝?huì)導(dǎo)致吉布斯現(xiàn)象太長會(huì)增加計(jì)算量。5. 二次開發(fā)與擴(kuò)展把這套代碼改成自己的工具5.1 封裝成可復(fù)用的函數(shù)原始代碼都是腳本形式變量直接暴露在工作區(qū)。如果你想反復(fù)使用某個(gè)功能比如生成任意參數(shù)的復(fù)指數(shù)序列可以把它封裝成函數(shù)function [n, x] gen_complex_exp(K, sigma, omega, n_start, n_end) % 生成復(fù)指數(shù)序列 x[n] K * exp((sigma j*omega)n) % 輸入 % K - 幅度系數(shù) % sigma - 實(shí)部衰減/增長因子 % omega - 虛部角頻率 % n_start - 起始索引 % n_end - 結(jié)束索引 % 輸出 % n - 時(shí)間軸向量 % x - 復(fù)指數(shù)序列 n n_start:n_end; x K * exp((sigma 1i*omega) * n); end調(diào)用時(shí)只需要一行[n, x] gen_complex_exp(1.5, 0.2, 3, 0, 15);。這樣你就不用每次改參數(shù)都去翻腳本了。5.2 用App Designer做個(gè)交互界面如果你想讓不懂Matlab的人也能玩這些實(shí)驗(yàn)可以用App Designer搭個(gè)簡單界面。拖幾個(gè)滑塊控制幅度、頻率、采樣間隔實(shí)時(shí)刷新波形圖。核心邏輯還是調(diào)用上面那些函數(shù)只是把輸入從命令行改成GUI控件。我一般會(huì)保留一個(gè)“導(dǎo)出數(shù)據(jù)”按鈕把當(dāng)前生成的信號(hào)和頻譜存成.mat文件方便后續(xù)分析。這個(gè)功能在寫論文或做項(xiàng)目時(shí)特別實(shí)用——你不需要重新跑一遍代碼直接加載數(shù)據(jù)就行。5.3 擴(kuò)展到實(shí)際信號(hào)處理這套實(shí)驗(yàn)用的是理想信號(hào)實(shí)際工程中你面對(duì)的是帶噪聲的音頻、傳感器數(shù)據(jù)或者通信信號(hào)。把實(shí)驗(yàn)里的采樣、頻譜分析、濾波這套流程套上去基本框架不用大改只需要注意幾點(diǎn)實(shí)際信號(hào)的頻率范圍未知先做頻譜預(yù)覽再定采樣率噪聲會(huì)讓頻譜變得毛糙需要加窗函數(shù)如hamming、blackman抑制泄漏濾波器設(shè)計(jì)要考慮實(shí)時(shí)性IIR濾波器階數(shù)低但相位非線性FIR濾波器線性相位但階數(shù)高5.4 驗(yàn)證方法改完代碼后怎么確認(rèn)沒跑偏我的習(xí)慣是拿一個(gè)已知答案的簡單信號(hào)做回歸測(cè)試。比如生成一個(gè)x[n] cos(0.1*pi*n)的序列理論上它的DFT應(yīng)該在k0.1*N處有一個(gè)峰值。跑一遍代碼看峰值位置對(duì)不對(duì)。如果對(duì)得上說明整個(gè)流程沒問題如果對(duì)不上就從采樣、加窗、FFT逐段排查。從那以后我每次改完信號(hào)處理代碼都強(qiáng)制走一遍這個(gè)“已知信號(hào)驗(yàn)證”流程。看起來多花了幾分鐘但省下了后面調(diào)試玄學(xué)問題的大把時(shí)間。希望這份實(shí)驗(yàn)記錄能幫你少走點(diǎn)彎路順利把信號(hào)分析與處理的實(shí)驗(yàn)跑通。本文還有配套的精品資源點(diǎn)擊獲取