信號時域采樣 頻譜分析(matlab)

基于matlab的時域信號采樣及頻譜分析一：主要設計方法與步驟：畫出連續(xù)時間信號x(t)=Ae-atsin(t)u(t)0因子A444.128，衰減因子a222.144，模擬角頻率 222.144；0x(tx(n)Aeantsin(nT)u(nn5001T 為采樣間隔，通過改變采樣頻率可改變T，畫出采樣頻率分別為200Hz，fs500Hz，1000Hz時的采樣序列波形；x(nx(t的幅頻曲線有無差別，如有差別說明原因；設系統(tǒng)單位抽樣響應為h(n)R5

(n)，求解當輸入為x(n)時的系統(tǒng)響應y(n)，畫出x(n)，h(n)，y(n)的時域波形及幅頻特性曲線，并利用結果驗證卷積定理的正確性（此內容將參數設置為A444.128，a222.144，W0222.144，fs1000；用FFT 對信號x(n)，h(n)，y(n)進行頻譜分析，觀察與4中結果有無差別；x(n)x1

(t)x1

(t)與原來的連續(xù)時間信號x(t)的時域波形，計算并記錄兩者最大誤差。二：詳細程序及仿真波形分析x(t及其200Hz/500Hz/1000Hzx(n)%繪制信號x(n)的幅度譜和相位譜clcclearallclosealln=0:50%定義序列的長度是50A=input('請入A的值A:')%設置信號的有關參數a=input('請入a的值a:')w0=input('w0T1=0.005T2=0.002T3=0.001T0=0.001x=A*exp(-a*n*T0).*sin(w0*n*T0)y1=A*exp(-a*n*T1).*sin(w0*n*T1)y2=A*exp(-a*n*T2).*sin(w0*n*T2)y3=A*exp(-a*n*T3).*sin(w0*n*T3)closeallsubplot(2,1,1)stem(n,x)%繪制x(n)的圖形gridontitle('離散時間信號')subplot(2,1,2)plot(n,x)gridontitle('連續(xù)時間信號')15010050

離散時間信號0-500 5 10 15 20 25 30 35 40 45 50連續(xù)時間信號150100500-500 5 10 15 20 25 30 35 40 45 50figure(2)subplot(3,1,1)stem(n,y1)gridontitle('200Hz理想采樣信號序列')subplot(3,1,2)stem(n,y2)gridontitle('500Hz連續(xù)時間信號')subplot(3,1,3)stem(n,y3)gridontitle('1000Hz連續(xù)時間信號')200100

理想采樣信號序列0-1000 5 10 15 20 25 30 35 40 45 50連續(xù)時間信號2000-2000 5 10 15 20 25 30 35 40 45 50連續(xù)時間信號2000-2000 5 10 15 20 25 30 35 40 45 50k=-25:25W=(pi/12.5)*kw=W/piY1=y1*exp(-j*pi/12.5).^(n'*k)figure(3)subplot(2,1,1)plot(w,abs(Y1))gridxlabel('w')ylabel('幅度')title('200Hz理想采樣信號序列的幅度譜')axis([-2201000])subplot(2,1,2)plot(w,angle(Y1))gridxlabel('w')ylabel('幅角')title('200Hz理想采樣信號序列的相位譜')1000度 幅

理想采樣信號序列的幅度譜0-2 -1.5

-0.5

0 0.5 w

1.5 2理想采樣信號序列的相位譜42角 0幅-2-4-2.5 -2 -1.5 -1 -0.5wY2=y2*(exp(-j*pi/12.5)).^(n'*k)figure(4)subplot(2,1,1)plot(w,abs(Y2))gridxlabel('w')ylabel('幅度')title('500Hz理想采樣信號序列的幅度譜')axis([-2201000])subplot(2,1,2)plot(w,angle(Y2))gridxlabel('w')ylabel('幅角')title('500Hz理想采樣信號序列的相位譜')

0 0.5 1 1.5 21000度 幅

理想采樣信號序列的幅度譜0-2 -1.5

-0.5

0 0.5 w

1.5 2理想采樣信號序列的相位譜42角 0幅-2-4-2.5 -2 -1.5 -1 -0.5wY3=y3*(exp(-j*pi/12.5)).^(n'*k)figure(5)subplot(2,1,1)plot(w,abs(Y3))gridxlabel('w')ylabel('幅度')title('1000Hz理想采樣信號序列的幅度譜')axis([-2201000])subplot(2,1,2)plot(w,angle(Y3))gridxlabel('w')ylabel('幅角')title('1000Hz理想采樣信號序列的相位譜')

0 0.5 1 1.5 21000度 幅

理想采樣信號序列的幅度譜0-2

-1

0 0.5 w

1.5 2理想采樣信號序列的相位譜42角 0幅-2-4-2.5 -2 -1.5 -1 -0.5w

0 0.5 1 1.5 21000Hz時沒有失真，500Hz時有橫線，產生失真，200Hz加長，失真加大。說明采樣頻率越大，失真越小。設系統(tǒng)單位抽樣響應h(n)R(nx(n)y(nx(n)，5h(n)，y(n)的時域波形及幅頻特性曲線，并利用結果驗證卷積定理的正確性（此內容將參A444.128a222.144，W0222.144fs1000。clcclearallclosen=1:50%定義序列的長度是50hb=zeros(1,50)%注意：matlab中數組下標從1開始hb(1)=1hb(2)=1hb(3)=1hb(4)=1hb(5)=1closeallsubplot(3,1,1)stem(hb)title('系統(tǒng)hb[n]')m=1:50%設定序列和長度值T=1%設定序列的采樣率A=1a=0.4T=1w0=2.0734x=A*exp(-a*m*T).*sin(w0*m*T)subplot(3,1,2)stem(x)title('y=conv(x,hb)subplot(3,1,3)stem(y)title('輸出信號y[n]')10.5

系統(tǒng)hb[n]00 5 10 15 20 25 30 35 40 45 50輸入序列x[n]10-10 5 10 15 20 25 30 35 40 45 50輸出信號y[n]10-10 10 20 30 40 50 60 70 80 90 100figure(2)subplot(3,1,1)plot(n,hb)gridontitle('矩形序列的時域波形')subplot(3,1,2)plot(x)gridontitle('輸入信號x[n]的時域波形')subplot(3,1,3)plot(y)gridontitle('輸出信號y[n]的時域波形')10.5

矩形序列的時域波形00 5 10 15 20 25 30 35 40 45 50輸入信號x[n]的時域波形10-10 5 10 15 20 25 30 35 40 45 50輸出信號y[n]的時域波形10-10 10 20 30 40 50 60 70 80 90 100分析：有數字信號處理中經常要進行卷積運算，conv可以用來計算兩個有限長序列卷積，該函數計算的兩個序列都是從n 0開始。用FFT 對信號x(n)，h(n)，y(n)進行譜分析，觀察與4中結果有無差別。clcclearallclosen=1:50hb=zeros(1,50)hb(1)=1hb(2)=1hb(3)=1hb(4)=1hb(5)=1closeallsubplot(3,1,1)m=1:50T=1A=1a=0.4T=1w0=2.0734x=A*exp(-a*m*T).*sin(w0*m*T)y=conv(x,hb)subplot(3,1,1)plot(n,abs(fft(hb)))title('h(n)FFT')subplot(3,1,2)plot(abs(fft(x)))title('x(n)的FFT')subplot(3,1,3)plot(abs(fft(y)))title('y(n)的FFT')64

h(n)的FFT200 5 10 15 20 25的FFT

30 35 40 45 502100 5 10 15 20 25的FFT

30 35 40 45 502100 10 20 30 40 50 60 70 80 90 100分析：matlabxDFT及其逆變換的函數分別為fft采用了混合算法，當N為質數時，采用的是原始的DFT算法。如果x用后計算出每列的NFFT。x(n)x1

(t)x1

(t)與原連續(xù)時間信號x(t)的時域波形，計算并記錄兩者最大誤差。%設置信號的有關參數clcclearallcloseA=input('pleaseinputtheA:')a=input('pleaseinputthea:')W0=input('pleaseinputtheW0:')fs=input('pleaseinputthefs:')n=0:49T=1/fst0=10/aDt=1/(5*a)t=0:Dt:t0xa=A*exp(-a*t).*sin(W0*t)K1=50k1=0:1:K1W1max=2*pi*500W1=W1max*k1/K1w1=W1/piXa=xa*exp(-j*t'*W1)x=A*exp(-a*n*T).*sin(W0*n*T)figure(1)subplot(4,1,1)plot(t*1000,xa)title('連續(xù)時間信號x(t)')axis([0t0*1000-50150])gridxlabel('t:毫秒')ylabel('x(t)')subplot(4,1,2)plot(w1,abs(Xa))title('連續(xù)時間信號頻譜Xa(w1)')axis([0100001200])subplot(4,1,3)stem(x)gridxlabel('n')ylabel('x(n)')title('采樣序列x(n)')axis([050-15160])x1=spline(n*T,x,t)gridxlabel('t:毫秒')ylabel('x(t)')subplot(4,1,4)plot(t*1000,x1)axis([0t0*10000200])title('x(n)gridxlabel('t:毫秒')ylabel('x1(t)')axis([045-20160])15010050x 0-500 5

連續(xù)時間信號x(t)15 20 25

35 40 4510005000

連續(xù)時間信號頻譜Xa(w1)150100t 50x 1501001 50x

0 100 200 300 400 500 600采樣序列x(n)0 5 10 15 20 25 30由0 5 10 15 20 25

700 800 900 100035 40 45 5030 35 40 45error=max(abs(x1-xa))k2=-25:25W2=(pi/12.5)*k2w2=W2/piX=x*(exp(-j*pi/12.5)).^(n'*kfigu