1、MATLAB Signal Processing Toolbox報告日期:2009/3/11指導老師:黃漢邦教授報告者:黃子豪OutlineTime response,ConvolutionFourier,Laplace,Z-transformDTFT and FFT IIR and FIR filtersSignal Processing Toolbox Example FDAtool Example波形產生波形產生Wihte noise-randn(t)Square function-square(t)Sinc函數-sinc(t)鋸齒波-sawtooth(t)步階響應與脈衝響應步階響應與脈
2、衝響應差分方程式y(n)+1.2y(n-1)+0.25y(n-2)=x(n),求出在n=0時的步階響應與脈衝響應filter指令格式 輸出信號(yn)=filter(輸入信號係數bm，輸出信號係數ak,輸入信號xn)NkMmmkmnxbknya00)()(步階響應與脈衝響應步階響應與脈衝響應clear;close all;t=-5:120;n0=0;%impulse starting pointa=1 1.2 0.25;%y(n)coefficientb=1;%x(n)coefficient%step signalx_step=(t=n0);%impulse signalx_imp=(t=n0
3、);%check responsestp_rep=filter(b,a,x_step);imp_rep=filter(b,a,x_imp);%plotfiguresubplot(211)stem(t,x_step);xlabel(time sequence(n);ylabel(input signal x);title(Step input signal)subplot(212)步階響應與脈衝響應步階響應與脈衝響應stem(t,stp_rep);line(t,zeros(1,length(t);xlabel(time sequence(n);ylabel(output signal y(n);
4、title(Step response)figuresubplot(211)stem(t,x_imp);xlabel(time sequence);ylabel(input signal x);title(Impulse input signal)subplot(212)stem(t,imp_rep);line(t,zeros(1,length(t);xlabel(time sequence(n);ylabel(output signal y(n);title(Impulse response)-2002040608010012000.51time sequence(n)input signa
5、l xStep input signal-20020406080100120-0.500.51time sequence(n)output signal y(n)Step response-2002040608010012000.51time sequenceinput signal xImpulse input signal-20020406080100120-2-1012time sequence(n)output signal y(n)Impulse response線性摺積運算線性摺積運算輸出性號y(n)可由將輸入訊號x(n)與系統脈衝響應T(n)作線性摺積而求得 函式:c=conv(
6、a,b)c的向量長度為length(a)+length(b)-1 範例:c=conv(1 1 1,1 1 1)c=1 2 3 2 1 0)()()(kknTkxny線性摺積運算線性摺積運算差分方程式y(n)+1.2y(n-1)+0.25y(n-2)=x(n),x(n)=0.5u(n),求出系統的響應%CHECK linear convolution AND full responseclear;close all;t=-5:20;n0=0;%impulse starting pointa=1 1.2 0.25;%y(n)coefficientb=1;%x(n)coefficient%step
7、signalx_step=(t=n0);x_in=(0.5.t).*x_step;%(0.5)n*u(n)%impulse signalx_imp=(t=n0);%impulse response stp_rep1=filter(b,a,x_imp);%convolutionstp_rep2=conv(x_in,stp_rep1);%輸入信號與脈衝響應的摺積%check response stp_rep3=filter(b,a,x_in);%使用filter指令求出系統響應與stp_rep2的結果比較%plotfiguresubplot(211)stem(t,x_in);xlabel(time
8、 sequence(n);ylabel(input signal x);title(Input signal)subplot(212)stem(t,stp_rep1);line(t,zeros(1,length(t);xlabel(time sequence(n);ylabel(impulse response y(n);title(Impulse response)線性摺積運算線性摺積運算線性摺積運算線性摺積運算figuresubplot(211)stem(t,stp_rep2(6:length(t)+5);line(t,zeros(1,length(t);xlabel(time seque
9、nce(n);ylabel(output signal y(n);title(Convolution)subplot(212)stem(t,stp_rep3);line(t,zeros(1,length(t);xlabel(time sequence(n);ylabel(output signal y(n);title(full response)-50510152000.51time sequence(n)input signal xInput signal-505101520-2-1012time sequence(n)impulse response y(n)Impulse respon
10、se-505101520-1-0.500.51time sequence(n)output signal y(n)Convolution-505101520-1-0.500.51time sequence(n)output signal y(n)full responseFourier transformVariable x(預設輸入變數符號),w(預設輸出變數符號)syms t w Inverse Fourier TransformVariable w(預設輸入變數符號),x(預設輸出變數符號)離散傅立葉轉換離散傅立葉轉換%sin signalclear;close all;w1=0.5*p
11、i;%analogue signal frequencyto=1/w1;%analogue signal periodts=0.5*0.5*to;%sampling periodtd=0:ts:4;%time sequencex_sin=sin(w1.*td);%fourier transform w=0:pi/50:2*pi;for i=1:length(w)for k=1:length(x_sin);x_ff(k)=x_sin(k)*exp(-j*w(i)*k);%傅立葉轉換公式 end x_f(i)=sum(x_ff(:);endnjwnjwenxex)(離散傅立葉轉換離散傅立葉轉換%p
12、lot sin figuresubplot(211)stem(td,x_sin);%plot the x-axis line(td,zeros(1,length(td);xlabel(time sequence n);ylabel(Sampling sin(nt);title(Sampling signal of Sin)subplot(212)plot(w/pi,abs(x_f);xlabel(pi unit);ylabel(Magnitude);line(w/pi,zeros(1,length(w/pi);title(Fourier transform of sin)00.511.522.
13、533.54-1-0.500.51time sequence nSampling sin(nt)Sampling signal of Sin00.20.40.60.811.21.41.61.82051015pi unitMagnitudeFourier transform of sin離散傅立葉轉換離散傅立葉轉換使用使用fft指令指令FFT是一種計算DFT的方法，可降低計算時間範例:clear%時間向量與訊號產生 t=(0:1/100:10-1/100);%Time vector x=sin(2*pi*15*t)+sin(2*pi*40*t);%Signal%FFT轉換 y=fft(x);%C
14、ompute FFT of x m=abs(y);%Magnitude p=unwrap(angle(y);%Phase(使相位連續)f=(0:length(y)-1)*99/length(y);%Frequency vector離散傅立葉轉換離散傅立葉轉換%畫出振幅圖 plot(f,m);title(Magnitude);%畫出相位圖 set(gca,XTick,15 40 60 85);figure;f=(0:length(y)-1)*99/length(y);plot(f,p*180/pi);title(Phase);set(gca,XTick,15 40 60 85);15406085
15、0100200300400500600Magnitude15406085-0.500.511.522.53x 104Phase逆離散傅立葉轉換逆離散傅立葉轉換-使用使用ifft指令指令%時間向量與訊號t=(0:1/255:1);x=sin(2*pi*120*t);figure;plot(t,x);title(x signal);%比較原始訊號x與重建訊號yy=real(ifft(fft(x);figure;plot(t,y);title(y signal);00.10.20.30.40.50.60.70.80.91-1-0.8-0.6-0.4-0.200.20.40.60.81y signal
16、00.10.20.30.40.50.60.70.80.91-1-0.8-0.6-0.4-0.200.20.40.60.81x signalLaplace TransformVariable t(預設輸入變數符號),s(預設輸出變數符號)Inverse Laplace TransformVariable s(預設輸入變數符號),t(預設輸出變數符號)Z TransformVariable n(預設輸入變數符號),z(預設輸出變數符號)範例:Syms z n f=n4;ztrans(f)ans=z*(z3+11*z2+11*z+1)/(z-1)5Inverse Z transformVariabl
17、e z(預設輸入變數符號),n(預設輸出變數符號)範例:syms z n f=2*z/(z-2)2;iztrans(f)ans=2n*n零點零點-極點分析極點分析)()1(.)2()1()1(.)2()1()(11zXznaazaaznbbzbbzYnanb%求解根b=2 3 4;a=1 3 3 1;zer=roots(b)pol=roots(a)k=b(1)/a(1)%化回多項式bb=k*poly(zer)aa=poly(pol)zplane(zer,pol)-1.5-1-0.500.511.5-1-0.500.513Real PartImaginary PartSignal Process
18、ing Toolbox Example 1-設計一帶通濾波器設計一帶通濾波器在command window下執行playshow filtdem1.三個頻率5,15,30Hz 的sine波Signal Processing Toolbox Example 12.設計一1020Hz帶通之IIR帶通濾波器Signal Processing Toolbox Example 13.5與30Hz sin波 被濾掉，剩下 15HzSignal Processing Toolbox Example 14.比較濾波前後的 頻率響應Signal Processing Toolbox Example 2-比較濾波
19、器差異比較濾波器差異在command window下執行playshow filtdem21.畫出理想的濾波器 頻率響應Signal Processing Toolbox Example 22.設計一YULEWALK 八階的 IIR filterSignal Processing Toolbox Example 23.設計一Butterworth and Chebyshev IIR filterSignal Processing Toolbox Example 24.頻率響應圖FDATool(Filter Design and Analysis Tool)MATLAB 命令視窗鍵入fdatoo
20、l指令fdatool1.選擇濾波器2.選擇濾波器設計方法3.指定濾波器階數4.輸入濾波器頻率規格5.輸入訊號衰減程度與漣波規格設計一最少階次帶通等漣波設計一最少階次帶通等漣波FIR濾波器濾波器取樣頻率(Fs)=2000Hz第一個阻帶結束頻率(Fstop1)=200Hz通帶開始的頻率(Fpass1)=300Hz通帶結束的頻率(Fpass2)=700Hz第二各阻帶開始的頻率(Fstop2)=800Hz第一各阻帶衰減大小(Astop1)=75dB通帶漣波大小(Apass)=0.1dB第二阻帶衰減大小(Astop2)=75dB帶通濾波，等漣波帶通濾波，等漣波FIR帶通帶通等漣波等漣波FIR濾波器濾波器
21、輸入帶通之頻率範圍輸入帶通之頻率範圍輸入振幅衰減的程度與漣波的振幅輸入振幅衰減的程度與漣波的振幅勾選勾選Minimum order按下Design Filter分析濾波器分析濾波器選擇上圖之分析工具列振幅圖相位圖將濾波器係數輸出到將濾波器係數輸出到workspaceFile-Export 將濾波器係數輸出成將濾波器係數輸出成C標頭檔標頭檔Targets Generate C Headerconst int BL=64;const int32_T B64=80690,-1072542,-2939985,5461255,6950698,-5912941,-1602054,-5104533,-111
22、44430,12287884,6241733,6077793,18813071,-23835745,-15696747,-5156699,-29284886,42299228,32893668,454314,43788106,-72822941,-65058054,12465656,-67152964,133647322,139847119,-51918454,129068678,-362675861,-580944589,712218408,712218408,-580944589,-362675861,129068678,-51918454,139847119,133647322,-671
23、52964,12465656,-65058054,-72822941,43788106,454314,32893668,42299228,-29284886,-5156699,-15696747,-23835745,18813071,6077793,6241733,12287884,-11144430,-5104533,-1602054,-5912941,6950698,5461255,-2939985,-1072542,80690;形成形成M-filefunction Hd=firbandstop%firbandstop Returns a discrete-time filter obje
24、ctFs=2000;%Sampling FrequencyFstop1=200;%First Stopband Frequency Fpass1=300;%First Passband FrequencyFpass2=700;%Second Passband FrequencyFstop2=800;%Second Stopband FrequencyDstop1=0.000177827941;%First Stopband AttenuationDpass =0.0057563991496;%Passband RippleDstop2=0.000177827941;%Second Stopba
25、nd Attenuationdens =20;%Density Factor%Calculate the order from the parameters using FIRPMORD.N,Fo,Ao,W=firpmord(Fstop1 Fpass1 Fpass2 Fstop2/(Fs/2),0 1 0,Dstop1 Dpass Dstop2);%Calculate the coefficients using the FIRPM function.b =firpm(N,Fo,Ao,W,dens);Hd=dfilt.dffir(b);set(Hd,Arithmetic,double);%Co
26、nvert the filter to the Direct-Form FIR structure.Hd=convert(Hd,dffir);set(Hd,Arithmetic,double);將濾波器係數輸出到將濾波器係數輸出到SPtool(Signal process tool)File-ExportSignal process tool命令視窗鍵入sptool訊號清單濾波器清單頻譜分析的清單使用FDAtool創建的帶通濾波器創立訊號來源創立訊號來源命令視窗輸入randn(state,0)x=randn(5000,1)SPToolFile-ImportSource from workspace變數x輸入為信號訊號變數x取樣頻率5000名稱Thanks for your attention
