計算機仿真技術(shù)(9)

1、1 概述 信號分類： 模擬信號 離散信號11 系統(tǒng)分類 1. 模擬信號處理系統(tǒng)(ASP)： （Analog Signal Processing） 模擬信號，即其幅值隨時間連續(xù)變化 模擬信號Xa (t)模擬信號處理器模擬信號ya (t)常用模擬信號類型： a.指數(shù)信號X(t)=Keat b.正弦信號X(t)=Ksin(t+) c.單位階躍信號 X(t)= 0 t=0 d.單位脈沖信號 0=t , t0 t=0 t001)(tX0)(t2. 數(shù)字信號處理系統(tǒng)(Digital Signal Processing) 由于DSP可以在通用的計算機上開發(fā)且可移植，這樣大大節(jié)約開發(fā)周期和成本，這是較ASP優(yōu)

2、越之處。但由于其運算速度低，特別對于強高頻信號，這是它的缺點。 12 DSP的研究內(nèi)容的研究內(nèi)容 (1) 信號分析時域分析時域分析示波器、最大值、最小值、均值、相關(guān)頻域分析頻域分析FFT分析、功率譜時頻域分析時頻域分析短時傅立葉變換（STFT）小波分析 復(fù)雜的信號可分解為簡單的信號的疊加，并以這些分量的組成來考慮信號的特征（幅值、相位）頻譜分析頻譜分析信號分析的基本任務(wù)：譜分析幅值譜 相位譜語音信號分析: 一維非平穩(wěn)時變圖像信號分析：二維自動測試控制系統(tǒng)：多信號 (2) 信號處理 對信號進行變換、加工、運算a. 濾波：把信號中最感興趣的部分提取出來b. 變換、調(diào)制、壓縮等 MATLAB交互式G

3、UI工具：sptool（信號觀察、濾波器設(shè)計以 及譜分析）13 信號的MATLAB表示 (1) 數(shù)字信號：按順序排列的一維信號采樣值向量序列如X4 3 7 -9 1行向量 XXX4 3 7 9 1 對于多通道采樣值用數(shù)組，每一列表示一個通道的信號1.4 波形的產(chǎn)生 (1) 首先產(chǎn)生一個時間向量（列向量） t0:0.001:1 ; (2) 產(chǎn)生一個50Hz和120Hz的正弦復(fù)合信號 y=sin(2*pi*50*t)+2*sin(2*pi*120*t); (3) 加上白噪聲信號 randn(state,0); yn=y+0.5*randn(size(t); plot(t(1:50);yn(1:50

4、)1.5 MATLAB特殊波形生成函數(shù) (1) 周期信號 sawtooth產(chǎn)生鋸齒波、三角波square產(chǎn)生方波【例】產(chǎn)生1.5秒50Hz的三角波，采樣頻率為10KHz，繪制0.2秒的波形圖 Fs=10000; T=0:1/Fs:1.5; X=sawtooth(2*pi*50*t); plot(T,X) axis( 0 , 0.2 , -1 , 1 ) Xmin,Xmax,Ymin,Ymax (2) 非周期信號 gauspuls產(chǎn)生高斯調(diào)制正弦脈沖 chirp產(chǎn)生線性快速掃描正弦函數(shù) pulstran脈沖串 sinc生成sinc函數(shù)dtcetj21)(sinxxxc)sin()(sin1.6

5、信號處理系統(tǒng)的描述（濾波器等） (1) 連續(xù)系統(tǒng)用微分方程描述(模擬濾波器) 系統(tǒng)函數(shù)(傳遞函數(shù)) 傳遞函數(shù)： 零極點模型：) 1()2() 1 () 1()2() 1 ()()()(11naaanbbbsUsYsHssssmmnn)()2()(1 ()()2()(1 ()()()(mpspspsnzszszsKsPsZsH狀態(tài)空間模型：(2) 離散系統(tǒng)模型(用于數(shù)字濾波器) a. 傳遞函數(shù) 注意：這種表示控制系統(tǒng)不同 b 零極點增益模型DuCXYBuAXXzzzzmnmaaanbbbzXzY) 1()2() 1 () 1()2() 1 ()()(11)()2()(1 ()()2()(1 ()

6、()()(mpzpzpznqzqzqzKzpzqzH c. 離散狀態(tài)空間 X(n+1)=AX(n)+Bu(n) Y(n)=CX(n)+Du(n) d 部分分式表達式 該式可用residuez得到znmKnmKKnpnrprzazbzzzn) 1()() 2 () 1 ()(1)() 1 (1) 1 ()()(11【例】 b=-4,8; a=1,6,8; r,p,k=residuez(b,a) r=-12 8 p=-4 -2 k=zzzzH21186184)(zzzH112184112)(e .二階因子級聯(lián)形式 Second-Order Sections(SOS) 式中，L為二階因子個數(shù) b01

7、 b11 b21 a01 a11 a21 k=1 b02 b12 b22 a02 a12 a22 k=2 SOS= b03 b13 b23 a03 a13 a23 k=3 b0L b1L b2L a0L a1L a1L k=LLkkkkkkkLkkzazaazbzbbHzzH122110221101)()( f. 格狀結(jié)構(gòu)(Lattice Structure)2. 信號的數(shù)字變換信號及系統(tǒng)分析的重要工具 2.1 連續(xù)信號的傅氏變換 (1) 周期信號的頻譜分析 原理：任何信號都可以分解成不同頻率成份的疊加。 數(shù)學工具：傅立葉級數(shù) 結(jié)論： a 周期信號的譜線是離散的 b 譜線只出現(xiàn)在基波頻率的整數(shù)

8、倍上 c 諧波分量幅值隨譜波次數(shù)的增加而減小(2) 非周期信號的頻譜分析 數(shù)學工具傅立葉變換 結(jié)論：頻譜是連續(xù)的*信號時域采樣dttfFetj)()()()()()(00 Xttx2.2 離散信號的傅氏變換 離散信號的表示：序列X(n) (1) 周期序列離散的傅氏級數(shù) 周期為N 令 離散傅氏級數(shù)系數(shù) DFS離散傅氏級數(shù))()(KNnXnXeknNjNnkXNnX210)(1)(10)()()(NnnknWkXnXDFSkX)(kX 可見，對于周期序列，離散傅氏級數(shù)系數(shù)也是周期序列。 系數(shù)周期序列是一個復(fù)數(shù)序列，應(yīng)在復(fù)頻域內(nèi)描述。 序列 序列傅氏級數(shù)系數(shù) 時序周期 在復(fù)頻域內(nèi)也是周期函數(shù) 周期N

9、 周期N 主值區(qū)間序列 主值區(qū)間序列 0,N-1 0,N-1WnknNkkXNkXIDFSnX10)(1)()(eWNjn2 )(nX)(kX 可見，對于周期序列，離散傅氏級數(shù)系數(shù)也是周期序列。 系數(shù)周期序列是一個復(fù)數(shù)序列，應(yīng)在復(fù)頻域內(nèi)描述。 序列 序列傅氏級數(shù)系數(shù) 時序周期 在復(fù)頻域內(nèi)也是周期函數(shù) 周期N 周期N 主值區(qū)間序列 主值區(qū)間序列 0,N-1 0,N-1(2) 有限長序列(非周期)離散傅立葉變換 0 n N-1 其余為0 有限長序列視為周期序列 的主值序列 離散傅立葉變換 離散傅立葉變換逆變換 X(k)是一個N長度序列，所對應(yīng)離散頻率=0-2 頻率間距為2 /N 0)()(nXnX

10、)(nXWnkNkkXNkXIDFSnX10)(1)()(10)()()(NnnkWkXnXDFSkXMATLAB函數(shù)：Xkdft(Xn，N)Xn=idft(Xk，N)2.3 快速傅立葉變換(FFT)注意到離散傅立葉變換enkNjNnnXkXnXDFT210)()()(WnkNnnX10)(eNjW2寫成矩陣形式：X(0) W00 W00 W00 x(0)X(0) = W00 W01 W(N-1)0 x(0)X(0) W00 W1 (N-1) W (N-1) (N-1) x(0)考察上述注意到：(1) Wnk的周期性：WWWnNkNknnk)()(eeNnkNjnkNj)(22(2) Wnk的

11、對稱性 基于上述W中，許多元素相同，使DFT 運算大大簡化。 eenkNjNnkNj2)2(2WWnkNnk2eeNNjnkNj2221ej方法：a 時間抽取基2FFT算法： X(n)：奇數(shù)項 N/2 FFT 偶數(shù)項 FFTb 頻率抽取基2FFT算法：前后對半分解FFT的MATLAB函數(shù)：fft( ) , ifft( )Y=fft(X)Y=fft(X,N)X序列： 單一信號：向量 多通道信號：矩陣每一列向量表示一通道信號， 并對每列向量求FFTYX的快速傅立葉變換N正整數(shù)，表示取N點的FFT若X向量長度N，則X補零；否則，截斷。Y一般為復(fù)序列：abs(Y)求幅值 angle(Y)求相位 頻率計

12、算：若X是在采樣頻率fs(Hz)下獲得，則Y序列兩點間的頻率為fs 則頻率軸：0:fs/2:fs Nyquist 頻率 0fs/2 和 fs/2fs 幅頻對稱性 用fft分析時，只考慮0 Nyquist 頻率 【例2.24】(1) 已知帶有測量噪聲信號 X(t)=sin(2f1t)+sin (2f2t) (t),其中，f150Hz，f2120Hz， (t)為均值為零的隨機信號，采 樣頻率1000Hz，采樣點數(shù)1024，用FFT繪制信號的頻譜圖，由該頻譜圖說明什么問題。%MATLAB PROGRAMFs=1000;N=1024;n=0:N-1;t=n/Fs;f1=50;f2=120;X=sin(

13、2*pi*f1*t)+sin(2*pi*f2*t)+2*randn(1,length(t)Y=fft(X,N)Mag=abs(Y);f=(0:length(Y)-1)*Fs/length(Y);plot(f(1:N/2),Mag(1:N/2) 可見，圖中在f50Hz，f120Hz處有兩個峰值，說明信號Y中含有兩個周期信號。 說明：在fft分析中，原信號序列X(n)長度Ndata和用于fft分析數(shù)據(jù)長度N，會對結(jié)果產(chǎn)生影響。一般取N為2 的冪Nn2，且取fft的N和Ndata相等避免因補零或截斷產(chǎn)生影響。2.4 隨機信號分析2.4.1 時域時域 自相關(guān)和互相關(guān)函數(shù)，隨機信號X(n)和Y(n) 自

14、相關(guān)： 互相關(guān)：)()()(mnXnXEmRx10)()(1NnNmnXnXNLim)()()(mnYnXEmRxy10)()(1NnNmnYnXNLim計算信號的自相關(guān)、互相關(guān)的MATLAB函數(shù)： c = xcorr(x) c = xcorr(x,y)自相關(guān)函數(shù)作用：識別隨機信號中是否含有相 同周期信號成份以及其頻率大小 互相關(guān)函數(shù)：工程測量中可找出兩個信號時 移大小。2.4.2 功率譜估計功率譜估計在頻域內(nèi)提取被噪聲淹沒的有用周期信號。自功率譜：自功率譜密度函數(shù)是自相關(guān)函數(shù)的傅立葉變換。互功率譜：由此，對于有限長隨機信號序列X(n) k=0,1,N-1 （求模） dfeRSfjxx2)(d

15、feRSfjxyxy2)()(21)(nXDFTSNfx方法：(1) 周期圖法(2) 分段平均周期圖法 a 將X(n)分成P個小段 b 對每段進行功率譜估計 c 求各段功率譜平均值作為整段功率譜值(3) 加窗分段平均周期圖法 對每段加非矩形窗 減小“頻率泄漏”，增加峰值窗度，提高估計分辨率。計算信號功率譜的MATLAB函數(shù)： psd自功率譜密度 csd互功率譜密度自功率譜自功率譜：Pxx=psd(x)Pxx=psd(x,Nfft,Fs,window,Noverlap)x信號序列NfftFFT長度Fs采樣頻率window窗函數(shù)，window=hanning(Nfft)Noverlap分段序列重疊

16、點數(shù)互功率譜互功率譜：Pxy=psd(x,y)Pxy=psd(x,y,Nfft,Fs,window,Noverlap)lW = baxcar(N)lW = hanning(N)lW = hamming(N)lW = bartlett(N)lW= blackman(N)lW= triang(N)lW= kaiser(N)lWchebwin(N,r)【例6.8】已知信號x(t)=sin(2f1t)+2sin(2f2t)+(t) 其中，f1=50Hz f2=120Hz (t)自噪聲， 信號長度1024， 采樣頻率Fs=1000Hz， 用函數(shù)psd繪制自功率譜密度估計曲線 。l%MATLAB PROG

17、RAM 6-8lclfl%Estimate PSD by Averaged Modified periodograml%Welchs Methodl%=lN=1024;lNfft=256;lFs=1000;ln=0:N-1;lt=n/Fs;lwindow=hanning(256);lnoverlap=128;ldflag=none;lxn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);lPxx=psd(xn,Nfft,Fs,window,noverlap,dflag);l%=l%Create frequency vectorl%=lf=(0:Nfft/

18、2)*Fs/Nfft;lsubplot(211)lplot(f,10*log10(Pxx)lxlabel(Frequency (Hz);lylabel(Power Spectrum (dB);ltitle(PSD-Welch Method)lgridl%Estimate PSD by Averaged Modified periodograml%Welchs Methodl%=lN=1024;lNfft=256;lFs=1000;ln=0:N-1;lt=n/Fs;lwindow=hanning(256);lnoverlap=128;ldflag=none;lxn=sin(2*pi*50*t)+

19、2*sin(2*pi*120*t)+randn(1,N);lPxx,f=psd(xn,Nfft,Fs,window,noverlap,dflag);lsubplot(212)lplot(f,10*log10(Pxx)lxlabel(Frequency (Hz);lylabel(Power Spectrum (dB);ltitle(PSD-Welch Method)lgridl 注意程序前半部分中頻率向量f的創(chuàng)建方法。這是與函數(shù)psd的輸出Pxx長度有關(guān)。 l 若x為實序列，l當Nfft為奇數(shù)時 f=(0:(Nfft+1)/2-1)/Nfft；l當Nfft為偶數(shù)時f=(0:Nfft/2)/Nff

20、t；1. 濾波概述 濾波是信號處理最基本、最重要、應(yīng)用最廣泛的技術(shù)。 濾波器分類： 模擬濾波器連續(xù)系統(tǒng) 數(shù)字濾波器離散系統(tǒng) 模擬 數(shù)字濾波器作用濾波器作用：(1) 對輸入信號的波形或頻譜進行加工處理， 以改善其特性。(2) 抑制甚至濾除某些頻率成份(噪聲，高頻信號)，提 取或恢復(fù)有用信號。數(shù)字濾波越來越多的用于解決工程實際問題： 數(shù)字音響 語言合成 通信 頻率合成 消除噪聲并從噪聲背景下提取有用信息 低通 lowpass按頻帶分類： 帶通 bandpass 高通 highpass 帶阻 bandstop 理想濾波器應(yīng)滿足的條件信號無失真?zhèn)鬏?信號無失真?zhèn)鬏敚盒盘柾ㄟ^濾波器后，其幅值與輸入信號成

21、比例關(guān)系，允許有時移，但不允許波形畸變。)()()()(jXjYKjH理想濾波器是物理不可實現(xiàn)濾波器；實際低通濾波器如下圖所示：濾波器設(shè)計性能指標：濾波器設(shè)計性能指標：通帶波紋 Rp(Passband ripple) dB阻帶衰減 Rs(Stopband attenuation) dB通帶邊界頻率 p 阻帶邊界頻率 s 過渡帶寬 (p s ) 實際濾波器是對理想濾波器幅值的逼進。通常濾波器幅頻特性|H(j)| 不是一個有理函數(shù)。因此用其平方幅頻響應(yīng)函數(shù)來代替)(22)(AjH 設(shè)計濾波器原型可化為其他三種濾波器。 MATLAB數(shù)字濾波器設(shè)計過程：2.1 原型濾波器設(shè)計函數(shù)原型濾波器設(shè)計函數(shù)(1

22、) 巴特沃斯濾波器(Butterworth) 模擬低通濾波器原型的平方幅頻響應(yīng)函數(shù)： 式中， c 為低通濾波器的截止頻率 N為濾波器的階數(shù)；N越大，過渡帶越小。 MATLAB設(shè)計函數(shù)： z,p,k=battap(N)cNAjH22211)()(【例】設(shè)計一個5階(或N=2,5,10,20)Butterworth低通模 擬原型濾波器，并繪制其平方幅頻響應(yīng)曲線。% n=0:0.01:2; N=5; z,p,k=buttap(N) b,a=zp2tf(z,p,k) H, =freqs(b,a,n) magH2=(abs(H).2; plot(,magHz)(2) Chebyshev 型濾波器 設(shè)計函

23、數(shù)： z,p,k=cheb1ap(N,Rp) Rp : 通帶波紋(3) Chebyshev 型濾波器 z,p,k=cheb2ap(N,Rs) Rs : 阻帶波紋(4) 橢圓濾波器 z,p,k=ellipap(N,Rp,Rs)(5) 貝塞爾濾波器 z,p,k=besselap(N) N: 階次cNCjHA222211)()(2.2 確定濾波器最小階次函數(shù)確定濾波器最小階次函數(shù) n, 1 =buttord(p , s ,Rp,Rs,s) n, 1 =cheb1ord(p , s ,Rp,Rs,s) n, 1 =cheb2ord(p , s ,Rp,Rs,s) s模擬濾波器2.3 頻率轉(zhuǎn)換函數(shù)頻率轉(zhuǎn)

24、換函數(shù)將低通模擬原型濾波器轉(zhuǎn)換為各類濾波器(低通，高通，帶通，帶阻) bt,at=la2lp(b,a, c ) (c :期望低通截止頻率) bt,at=la2hp(b,a, c ) (c :期望高通截止頻率) bt,at=la2bp(b,a, o ,Bw ) bt,at=la2bs(b,a, o ,Bw ) o ：中心頻率 Bw：帶寬， 下 上 2.4 模擬濾波器完全設(shè)計函數(shù)模擬濾波器完全設(shè)計函數(shù) (1) 巴特沃斯模擬濾波器 b,a=butter(n, o ,ftype,s) n階次 o截止頻率 high ：高通 ftype類型： stop ：帶阻 缺?。旱屯ɑ驇?s模擬【例】設(shè)計一個But

25、terworth模擬帶通濾波器，設(shè)計指標：通帶頻率10002000Hz，兩側(cè)過渡帶寬500Hz， 通帶波紋：1dB，阻帶衰減大于100dB。%MATLAB PROGRAM 3-11%Design a Butterworth Analog bandpass filter%Desired performents of the filterwp=1000 2000*2*pi;ws=500 2500*2*pi;Rp=1;Rs=100;%Compute Oder and Cutoff frequencyN,Wn=buttord(wp,ws,Rp,Rs,s);NFc=Wn/(2*pi)%Compute t

26、he filterb,a=butter(N,Wn,s);%Outputw=linspace(1,3000,1000)*2*pi;H=freqs(b,a,w);magH=abs(H);phaH=unwrap(angle(H);plot(w/(2*pi),20*log10(magH);xlabel(Frequency (Hz);ylabel(Magnidute (dB);title(Butterworth Analog bandpass filter)grid on(2) 其他直接設(shè)計函數(shù)其他直接設(shè)計函數(shù) b,a=cheby1(n,Rp, n ,ftype,s) b,a=cheby2(n,Rs,

27、n ,ftype,s) b,a=ellip(n,Rp,Rs, n ,s) b,a=besself(n, n ,ftype)2.3 數(shù)字濾波器設(shè)計數(shù)字濾波器設(shè)計(1) 類型a 無限沖激響應(yīng)數(shù)字濾波器(IIR) Infinite Impulse Response Digital Filter010)(1)()()(nnkNkkrMrrzzazbnhzXzYzHb 有限沖激響應(yīng)數(shù)字濾波器(FIR) Finite Impulse Response Digital Filter h(n)濾波器脈沖響應(yīng)有限項 按頻率特性： 低通 高通 帶通 帶阻zbkNnkzHzXzY10)()()(2) IIR濾波器的經(jīng)典設(shè)計法 借助于模擬濾波器設(shè)計 離散化 IIR濾波器 離散化函數(shù)離散化函數(shù)： bz,az=impinvar(b，a，F(xiàn)s)脈沖響應(yīng)不變法 b，a模擬濾波器 Fs采樣頻率 bz,az=bilinear(b，a，F(xiàn)s)雙線性變換 IIR濾波器完全設(shè)計函數(shù)濾波器完全設(shè)計函數(shù)： b,a=butter(n，n ,ftype) 注意： 后面無s n : 標準化數(shù)字頻率01.0【例】設(shè)計一個butterworth高通數(shù)字濾波器，通帶邊界頻率為300Hz，阻帶邊界頻率200Hz，通帶波紋1dB，阻帶衰減大于20d