西南交通大學(xué)信號處理期末作業(yè)

1、1、考慮兩個諧波信號x(t)和y(t),其中Mf) = Acos(wct +, y(t) = Bcos(vv/)式中4和久為正的常數(shù),。為均勻分布的隨機變量，其概率密度函數(shù)為/) =0,其他而8是一個具有零均值和單位方差的標準高斯隨機變量，即其分布函數(shù)為fB(b)=4=exp(-/72 / 2), yc b 8Q2兀(1)求W)的均值“)、方差b、自相關(guān)函數(shù)凡和自協(xié)方差函數(shù)q(r)0(2)若。與6為相互統(tǒng)計獨立的隨機變量，求.*)和.*)的互相關(guān)函數(shù)弋/)與互協(xié)方差函 數(shù) J ”)。解：x(f)的均值，為： TOC o 1-5 h z 1 r”1=0o儲=E(j(/) = E(Acos(vvf

2、z + 0)= Acos(wet + 族卜/0 =+ ) 2tt 為2兀方差。為:.,.,AAAa； = E(x2(t) = E(A2 cosz(vv.r + 姆)=E(l + cos(2vv./ + 20)= + E(cos(2卬/ + 2。)= 自相關(guān)函數(shù)6(r)為：/?)A cos(wet+wcr + )= A:E(cos(h;7 +(/) cos(wct+wc T +(/)A 二A2AA=,(cos(2vv;/4-2vvrr + 2 姆 + cos(u;.r)= cos(H；,r) + (cos(2vvrZ+2iv,r + 2)= cos(vvr) 自協(xié)方差函數(shù)耳(r)為：c*) =

3、R) = cos(ivr)y(f)的均值為：uy (t)=嗎(y(7)= EB(Bcos(wct)= EB)cos(vv/) = 0，所以 E(B)=O 由互相關(guān)函數(shù)的定義可知：R. (r) = E(Acos(h;/ + Q)Bcos(町f- wcr) 由題意知道。與8為相互統(tǒng)計獨立的隨機變量，所以有R、(t)= *(4cos(w/ + e)X8cos(*f-+T) = 4cos(w/ +e)E(8)cos(w/ 卬/) =A x 0 + 0 x cos(vv/ vv. r) = 0 互協(xié)方差函數(shù)c,(D%(r)=q(r)=。.接收信號由下式給出：y = Acos(qi + d) + q,i

4、= L2,.,N ,式中弓%。1)即是零均 值和單位方差的高斯噪聲，為載波角頻率，而8是未知的相位。假設(shè)如牡，。相互獨立, 求未知相位的最大似然估計。解：由于如在，。相互獨立，所以)，.以也相互獨立并且服從高斯分布，可以得到 曲.以與8的聯(lián)合戒率密度函數(shù)分布由此，可以得到似然函數(shù)N N_L = -ln(2r)- -22yr- 4 cos(g/ + 0)22 i-i該似然函數(shù)對e求偏導(dǎo)，并令該偏導(dǎo)函數(shù)為零，即可得到如下公式: 丸 n=人 cos(例 i + e)sin(i+e)= O該函數(shù)為非線性方程不容易求解,Z A cos(、i + Oml) sin(卬 + 6ML) = Z 凹 sm(、i

5、 + Oml )因此，最大似然估計6M.為上述函數(shù)的零點值。 則若忽略雙倍頻率2 ,則可簡化到如下式子:NAZysin(Qj+e)= o根據(jù)三角公式分解得到如下式子:A jVsm6ML Z凹 cos(如)=-cos% g)：由此，可以得到如下公式N- Zxsm&itan &ml = Z y cos(Q.i)1-1所以相位的最大似然估計如下：NAZ E sm(ei)般=arctan(Z cosgi).離散時間的二階AR過程由差分方程x() = 4M-1) +3(-2) +叭)描述，式中貝) 是一零均值、方差為。:的白噪聲。證明式)的功率譜為bp(f)=；1 + a； + 6/； 24 (1 -

6、% ) cos(2兀/) - 2a2 cos(41 f )證明：由AR過程的功率譜公式知(/)=-_廿-其中1 -azej4JT, =(-aYeiZ!tf azeJ4irf)(1 -aleJ2ir/ - a?/)=+ a； + a； _ a二(/ 4- e41Tf ) + alazeiZftf + ala2e)2JTf al(ej2r/ + eiZxf )=1 + a：-2a1(l-a2)cos(2/)-2az cos(4f)將其帶入第一個公式可得：a2?=；21 + a； + a； - 2% (1 - % ) cos(2兀/) - 2a2 cos(41 f )信 號 的 函 數(shù) 表 達 式

7、為：X(/) = sin(2100/) +1.5sm(2300/) + A(t)sin(2200/) + dn(t) + n(t),其中，4(f)為一隨時間變 化的隨機過程，疝為經(jīng)過390-410HZ帶通濾波器后的高斯白噪聲，(，)為高斯白噪聲， 采樣頻率為1kHz,采樣時間為2.048s。分別利用周期圖譜、ARMA、Buig最大燧方法估計 信號功率譜，其中ARMA方法需要討論定階的問題。解：由題意知采樣點數(shù)一共為：1000 x2.048=2048個數(shù)據(jù)點。A(f)為一隨時間變化的隨 機過程，由于隨機過程有很多類型，如維納過程、正態(tài)隨機過程，本文采用了均值為0,方 差為1的正態(tài)隨機過程來作為演

8、示，來代替A。)，高斯白噪聲采用強度為2的高斯白噪聲 代替(/)，其帶通濾波后為力小)。其中濾波器采用的是契比雪夫數(shù)字濾波器。 可得到x(t)如下圖所示：1、周期圖法matlab中的周期圖功率譜法原理是通過計算采樣信號的FFT,獲得離散點的幅度，再根 據(jù)幅度與功率之間的關(guān)系，轉(zhuǎn)換為離散點的功率，再通過坐標變換將離散點的功率圖轉(zhuǎn)換為 連續(xù)功率譜密度。S/cp/:計算采樣信號x(ii)的DFT,使用FFT方法來計算。如果此處將及頻率處的幅度對 稱到物理實際頻率，得到的就是單邊譜，否則就是雙邊譜根據(jù)正余弦信號功率與幅度的關(guān)系以及直流功率與幅度的關(guān)系，將幅度轉(zhuǎn)換為離 散功率譜。對橫縱坐標進行轉(zhuǎn)換，橫坐

9、標乘以頻率分辨率轉(zhuǎn)換為實際連續(xù)物理頻率，縱坐標 除以頻率分辨率轉(zhuǎn)換為功率譜密度。調(diào)用 MATLAB 中自帶的 matlab 中Pxx.f=periodog】am(x,window.nfitfs)函數(shù)可得計算結(jié) 果如下：周期圖法求功率譜 TOC o 1-5 h z 21.5 -玄0050100150200250300350400450500f/fs5.附件中表sheet1為某地2008年4月28 口凌晨12點至2008年5月4 口凌晨12點的 電力系統(tǒng)負荷數(shù)據(jù)，采樣時間間隔為1小時，利用Kalman方法預(yù)測該地5月5 口的電力系 統(tǒng)負荷，并給出預(yù)測誤差(5月5 口的實際負荷數(shù)據(jù)如表sheet2)

10、。解：卡爾曼濾波是以最小均方誤差作為估計的最佳準則，來尋求一套遞推估計的算法, 其基本思想是：采用信號與噪聲的狀態(tài)空間模型，利用前一時刻地估計值和現(xiàn)在時刻的觀測 值來更新對狀態(tài)變量的估計，求得出現(xiàn)時刻的估計值。它適合于實時處理和計算機運算。現(xiàn)設(shè)線性時變系統(tǒng)的離散狀態(tài)防城和觀測方程為：X(k)=F(k,k-1 )X(k-1 )+T(k,k-1 )U(k-1)Y(k) = H(k)X(k)+N(k)其中：X(k)和Y(k)分別是k時刻的狀態(tài)矢量和觀測矢量，F(xiàn)(k.k-1)為狀態(tài)轉(zhuǎn)移矩陣，U(k)為k 時刻動態(tài)噪聲，T(k,k-1)為系統(tǒng)控制矩陣，H(k)為k時刻觀測矩陣，N(k)為k時刻觀測噪聲。

11、 卡爾曼濾波的算法流程為： 1、預(yù)估計X(k) =F(kl)X(k-l)2、計算預(yù)估計協(xié)方差矩陣e(k)=F(k,k-l)xC(k)xF(k,k-D+T(k,k-l)xQ(k)xT(k,k-l)，Q(k)=U(k)xU(kY3、計算卡爾曼增益矩陣K(kC(k) xH(k)rxH(k)x C(k) xH(k)f+R(k)-R(k尸N(k)xN(k)f4、更新估計X(k)=X(k) +K(k)xY(k)-H(k)x X(k)5、計算更新后估計協(xié)方差矩陣C(k) = LK(k)xH(k)x C(k) xi-K(k)xH(k)T+K(k)xR(k)xK(k)，X(k+l)=X(k)C(k+l)=C(k

12、)6、重更以上步驟最終可以獲得如下結(jié)果：使用Kalrmm對電力系統(tǒng)負荷數(shù)據(jù)進行頸測35003000250020001500100050020406080100120140160180時間點數(shù)350030002500使用Kalman對電力系統(tǒng)負荷數(shù)據(jù)進行預(yù)測300200 1000 100預(yù)測值與真實值之誤差168170172174176178180182184186188190192時間點數(shù)本題將表中的作為觀測數(shù)據(jù)，圖中橫坐標為1表示2008.4.28 1時刻數(shù)據(jù),2表示2008.4.28 口 2時刻的數(shù)據(jù)，一次類推，168表示200855 口 1時刻的數(shù)據(jù)。從表中可以看出預(yù)測誤差 的最大值為3

13、00。預(yù)測誤差的大小與代碼中的R、Q值得設(shè)置有關(guān)。Q越大預(yù)測誤差越小， 但是同時也表明系統(tǒng)內(nèi)的噪聲很大。本題中取得Q、R值均為高斯分布的協(xié)方差。代碼見附錄。6,設(shè)某變壓器內(nèi)部短路后，故障電流信號分解得到卜.式:y(0= 20et/x + 20sin (ent + 60 ) + 12sin(2cnt + 45 ) + lOsin(3(nt + 30 ) + 6sin(4(ot + 22.5 ) + 5sin(5wt + 36 )式中，fn = 27ll0, T 30ms,o = 5OH/分別利用小波變換、短時傅里葉變換和維格納威利分布分析故障電流信號的時頻特性。解：（1）小波變換：連續(xù)小波變換的

14、定義：。啊(“,s) =匚/(0 爐(彳)力計算連續(xù)時間小波變換的4個步驟:選取一個小波，然后將其和待分析信號從起點開始的一部分進行相乘積分。 計算相關(guān)系數(shù)C。將小波向右移，重復(fù)U和2的步驟直到分析完整個信號。將小波進行尺度伸縮后再重及1, 2, 3步驟，直至完成所有尺度的分析。（2）短時傅里葉變換短時傅里葉變換定義如下:STF/ (,) = /(/), g“=匚 /(f)gQ - )產(chǎn)力STFTf （,9=/OgQ37（助（出一9*吐文/口J24j（3）維格納威利分布變換 維格納威利分布定義如下：WDt（/.Q）= （： / + 弁r在MATLAB中沒有維格納威利分布變換的相關(guān)函數(shù)，需要安裝

15、一個MATLAB版本的時頻 分析工具箱。調(diào)用里面的函數(shù)即可。小波變換和短時傅里葉變換MATLAB均自帶了相關(guān)的 函數(shù)。程序見附錄。代碼運行結(jié)果結(jié)果如下：小波時頻圖501001502002500.10.20,3040.50.60.70 80.9時間t/S70605040302010妞時停單葉交換結(jié)果50000.10.20.30.40.50.60.70.80.91時間t/s7.假定一電力系統(tǒng)諧波與間諧波信號的函數(shù)表達式如下：y () = 0.001 cos(2 乃 x 10 + 疚)+ cos(2 乃 x 50 + 4)+ 0.1 cos(2;r x 150 + 內(nèi))+ 0.002 cos (2乃

16、 x 50 + 四)+()其中，采樣頻率為1024Hz，相位-內(nèi)為獨立的均勻分布-,+;4)為一噪聲信號，信噪比取為20dB ,分別采用三種現(xiàn)代信號處理方法進行諧波與間諧波頻率提取與譜估計。解:本題目采用的頻率提取的三種方法為小波變換、短時傅里葉變換和維格納威利分布。采用周期圖法、MUSIC法、Burg法進行譜估計。確定出諧波的頻率為50Hz和150Hz。原始侑號200.20,4060.811.21.41.61.8時向（U秒） 小波時頻圖00.20.40.60.811.21.41.61.8時間t/S 500 400 300 200 100 0Wignor-Villo timo-froquonc

17、y distribution500450400350z 300三2502001501OO500.20.40.60.811.2位同Vs1.41.61.80.5周期圖法求功率譜4 3 2 1 o o o oqp、w-a50100 150 200 250 300 350 400 450 500 550 600f/HzMUSIC方法o o o o4 2 2- (gp)與-40050100150 200 250 300350 400 450500頻率(f/Hz)Burg法譜估計 0.0150.010.0050 050 100 150 200 2500060附錄代碼：第四題：clc;clear;fs=10

18、00;%采樣頻率T=2.048;%采樣時間t=0:l/fs:T;A = normmd(0,l JJength(t);%方差為 1,均值為。的高斯分布N=wgn(l ,length(t),2);%強度為2的高斯白噪 聲300 350 400 450 500 550 600f/fsDn=bandp(N,390,410,200,450,0.1,30.fs);figure。)；subplot(211);plot(t,N);title(,原始高斯白噪聲工subplot(212);plot(t.Dn);titleC帶通濾波后高斯白噪聲)；Sig=siii(2*pi* 100. *t)+1.5 *sin(2

19、 *pi*300.*t)+A.*siii(2 *pi*200.*t)+Dn+N;figure。)；plot(t,Sig);titlef原始輸入信號)axis(O 2.1 -7 7);%周期圖譜Pxx,f=periodogiam(Sig,lengtli(t),fs);% 周 期圖法figuie(3);plot(f,Pxx);出le(，周期圖法求功率譜工xlabelCf/Hz); ylabelC 功率/db);% ARMA譜估計z=iddata(Sig);%#信號轉(zhuǎn)化為matlab接受的 格式RecoidAIC=;fbrp=l:20%自回歸對此PACF,給定滯后長 度上限p和qfor q=l:20

20、%移動平均對應(yīng)ACF m=ainiax(z(l :lengtli(t).p,q);AIC = aic(m); %armax(p,q)選擇對 應(yīng)FPE最小，AIC值最小模型RecordAIC=RecordAIC;p q AIC;end endfbr k=l:size(RecoidAIC,l)ifRecoidAIC(k,3)=niui(RecoidAIC(：3) %選 擇AIC最小模型pa_AIC=RecordAIC(k, 1);qa_AIC=RecordAIC(k,2); break;end endmAIC=aimax(z( 1 :length(t). pa_AIC.qa_AIC );Pxx2,

21、f2=freqz(niAIC.C4iiAIC.a,fs);P2=(abs(Pxx2).* 1).A2;P2tol=10*logl0(P2);figuie(4);plot(f2/pi*fs/2,P2tol); title(ARMA 法(AIC 準 則)Xabul(Hz);ylabulC 振幅/dB);plot(RecordAIC(: ,3);ylabel(AIC(p,q),);% burg法計算Pxx,F = pburg(Sig,60,length(t),fs);%burg 法 figure(6);plot(EPxx);titlefBuig法譜估計工xlabel(,f/fs,);%X 軸坐標名稱

22、ylabelC功率譜/dB);%Y軸坐標名稱 %fiinction y=bandp(x,f Lf3,fsl,Ah,ip,rs.Fs)%帶通濾波%使用注意事項：通帶或阻帶的截止頻率與 采樣率的選取范圍是不能超過采樣率的一 半%即,fl,B,fsl,fsh,的值小于 Fs/2%x:需要帶通濾波的序列% f 1 ：通帶左邊界% f 3：通帶右邊界%fsl：衰減截止左邊界% fsh：衰變截止右邊界%ip：邊帶區(qū)衰減DB數(shù)設(shè)置%rs：截止區(qū)衰減DB數(shù)設(shè)置%FS：序列x的采樣頻率% fl=3OO;f3=5OO;%通帶截止頻率上下限% %=200；fsh=600;%阻帶截止頻率上下限% ip=0.1;rs=

23、30;%通帶邊衰減DB值和阻帶 邊衰減DB值%Fs=2000;% 采樣率 %wpl=2*pi*fl/Fs;wp3=2*pi*f3/Fs;wsl=2*pi*fsLTs;wsh=2*pi*fsh/Fs;wp=wpl wp3;ws=wsl wsh;%設(shè)計切比雪夫濾波器；n,wn=cheb 1 ord(ws/pi,wp/pi,rps);bz 1 ,az 1 =cheby 1 (n,rp,wp/pi);%查看設(shè)計濾波器的曲線hjv=fieqz(bz 1 ,az 1,256,Fs);h=20*log 10(abs(h);y=filter(bzl.azLx);end第5題%本題目需要提醒一點：給的數(shù)據(jù)為觀測

24、數(shù) 據(jù)Z而不是X clc;clear;xl=xlsiead(/負荷數(shù)據(jù).xlsYshee); xl=xl(:,2);x2=xlsieadC./負荷數(shù)據(jù).xls?sheet2);x2=x2(:,2);x=xl;x2;Nl=length(xl);N=length(x);A=l；B=0;H=l；w=nomund(0,1000,1 ,N);% 這里隨便取值 v=noimrnd(0.1000,1 ,N);P(1尸16;%隨便取值Z=x;X(l)=24;%隨便取值R=cov(v);Q=cov(w);fbr 1=2 :Ntempx=A*X(i-1 );%+B *u ； TempP=A*P(i-1 )* A+

25、Q;K(i)=TempP*H,* 1 /(H*TempP*Hr+R); X(i)=X(i-1 )+K(i)*(Z(i)-tempx);P(i)=( 1 -K(i)*H)*TempP;endt=l:length(Z);figure;plot(t,Z,b;t,X(t),T);titleC使用 Kahuan對電力 系統(tǒng)負荷數(shù)據(jù)進行預(yù)測工xlabelf時間點數(shù));ylabe(電力系統(tǒng)負荷 *);axis tight;legend(，負荷真實值？Kalman 預(yù) 測值工 figure;subplot(2,l,l);t=length(x 1): length(x);使用 Kalman 對電 力系統(tǒng)負荷數(shù)據(jù)

26、進行預(yù)測力xlabelC時間點數(shù) );ylabel(電力系統(tǒng)負荷力axis tight;legend(,負 荷真實值？Kalman預(yù)測值)； set(gca,XTick,length(xl):2:length(x); subplot(2,l,2);enoi-Z-X*;plot(t,enor(t);titleC預(yù)測值與真實值之誤差 力xlabelC時間點數(shù)) set(gca/XTicklength(x l):2:length(x);ylabel(,5月5 口預(yù)測值與真實值誤差axis tight;第六題：%小波變換clc;clear;close all;f=50;%信號頻率ouniiga=2*p

27、i*f;N_sample=2048;%總采樣點數(shù)Fs=1000;%采樣頻率 t=O:l/Fs:l;Tao=0.03;A=l;%信號幅度X=20*exp(-t/Tao)+20*sm(ouniiga*t+pi/3)+12*si n(2 *ounuga *t+pi/4)+10*siii(3 *ounuga*t+piz6 )+6*sin(4*oumiga *t+pi/8)+5 * sm(5 *ouiniga *t +pi/5); % 信號函數(shù)表達式figure;plot(t,x);titl原始信號工xlabelfWfu t/sTontSize14);ylabel(幅值；FontSiz-M);%原信號函

28、數(shù)wavenaine=,cmor3 -3 totalscal=256;Fc=centfrq(waveiiaine); %小波中心頻率 c=2*Fc*totalscal;scals=c./(l :totalscal);fscal2fiq(scals,wavename,l/Ts); % 將尺度 轉(zhuǎn)換為頻率coefs=cwl(x,scals,wavename); % 求連續(xù)小 波系數(shù)figure;unagesc(tXabs(coefs);colorbar;xlabelfWfu t/sTontSize14);ylabel(頻率 f HzTontSize 14);title (,小波時頻圖；FontSi

29、ze, 16);axis(0 1 0 300);%短時傅里葉變換S,F,T,P=spectiogiam(x,256,250,256,Fs);figure;surf(TE10*logl0(P)/edgecolor7none,)； axis tight;view(0,90);xlabelfRt 間/s); ylabeiq頻率/Hz);titleC短時傅里葉變換結(jié)果力% Wigner-Ville time-frequency distnbution. X=hilbert(xf);scals=c./(l :totalscal);fscal2fiq(scals,wavenameJ/Ts); % 將尺度

30、轉(zhuǎn)換為頻率coefs=cwt(x 1,seals,wavename); % 求連續(xù)小波系數(shù)figure;unagesc(n,f,abs(coefs);colorbar;xlabelC時間 t/s,FontSize, 14);ylabel(瀕率 f;HzTontSize 14);titleC 小波時頻圖；FontSize： 16);figure;contour(t/Fs,#Fs,abs(tft);xlabelC時間 t/sr);ylabelC頻率 fHz*);title(,Wigner-Villetiine-fiequencydistribution); axis(0 1 0 300) % 第七題： clc;clear; close all;%參數(shù)設(shè)置Fs= 1024; %采樣頻率 n = 0:LFs201;%采樣時間 N = length(n); % 采樣點Wl=0.001 *cos(2 *pi*n*l O+unifhid(-pi,pi)+c os(2*pi*50*n+unifrnd(-pi,pi)+0.1 *cos(2*pi*