1、考慮兩個(gè)諧波信號(hào) 和 ，其中 , 式中 和()xtyt()cos()xtAwt(cos()ytBwtA為正的常數(shù), 為均勻分布的隨機(jī)變量，其概率密度函數(shù)為cw,1,02()f其 他而 是一個(gè)具有零均值和單位方差的標(biāo)準(zhǔn)高斯隨機(jī)變量，即其分布函數(shù)為B 21()exp(/),Bfbb（1）求 的均值 、方差 、自相關(guān)函數(shù) 和自協(xié)方差函數(shù) 。()xt()xut2xtxR()xc（2）若 與 為相互統(tǒng)計(jì)獨(dú)立的隨機(jī)變量，求 和 的互相關(guān)函數(shù) 與互協(xié)方差 ()t yR函數(shù) 。xyc解：（1）的均值 為：()t()xut 220 011(cos)cos()sin() 2cEAwtAwtdAwt方差 為：2()xt2 22(cs)(1cs)(os2)c c cAtEtEt自相關(guān)函數(shù) 為：)xR22 2 2()os(o(+)(o)s(+)c+2)cscsc2cos()xc cEAwttwAwttww 自協(xié)方差函數(shù) 為：()x 2()os()xxcR（2） 的均值為：()yt,所以()()cos()()0yBBcutEtwtEBt()=0EB由互相關(guān)函數(shù)的定義可知： osxyccRAw由題意知道 與 為相互統(tǒng)計(jì)獨(dú)立的隨機(jī)變量，所以有()cos()(cos)()(s)00xy ccccRAwtttt 互協(xié)方差函數(shù) xy ()0xyxyR2.接收信號(hào)由下式給出： ，式中 即 是零cos,12,.i iAiN(0,1)ii均值和單位方差的高斯噪聲， 為載波角頻率，而 是未知的相位。假設(shè) 相互2.N獨(dú)立，求未知相位的最大似然估計(jì) 。ML解：由于 相互獨(dú)立，所以 也相互獨(dú)立并且服從高斯分布，可以得到12,.N1,.Ny與 的聯(lián)合概率密度函數(shù)分布1,.Ny 21cos()1(,.|)(2)iyAiNNfye由此，可以得到似然函數(shù)21ln(2)cos()NiiLyAi該似然函數(shù)對(duì) 求偏導(dǎo)，并令該偏導(dǎo)函數(shù)為零，即可得到如下公式：1cos()sin()0Ni ciy 因此，最大似然估計(jì) 為上述函數(shù)的零點(diǎn)值。ML則 1 1cos()sin()sin()N NMLMLMLc ci iAi y 該函數(shù)為非線性方程，不容易求解，若忽略雙倍頻率 ，則可簡(jiǎn)化到如下式子：2c1si()0Nciy根據(jù)三角公式分解得到如下式子： 1 1sinos()ssin()NNMLicMLciiy由此，可以得到如下公式 1sitanco()NiMLiy所以相位的最大似然估計(jì)如下： 1sin()arct(oNciMLiy3.離散時(shí)間的二階 AR 過(guò)程由差分方程 描述，式中12()()(xnaxnw是一零均值、方差為 的白噪聲。證明 的功率譜為()wn2w22112()()cos()cos(4)wxPfaafaf證明：由 AR 過(guò)程的功率譜公式知 2241()xjfjfPfae其中 22424241 114 2221122()() (coscosjfjf jfjfjfjfjffjfjfjfjfae aeeaee將其帶入第一個(gè)公式可得： 21122()()cs()cos(4)xPfaafaf4、信號(hào)的函數(shù)表達(dá)式為：，其中， 為一隨時(shí)間變 sin(210).5sin(230)sin(20)xtttAtdntAt化的隨機(jī)過(guò)程， 為經(jīng)過(guò) 390-410Hz 帶通濾波器后的高斯白噪聲， 為高斯白噪聲，d n采樣頻率為 1kHz，采樣時(shí)間為 2.048s。分別利用周期圖譜、 ARMA、Burg 最大熵方法估計(jì)信號(hào)功率譜，其中 ARMA 方法需要討論定階的問(wèn)題。解：由題意知采樣點(diǎn)數(shù)一共為：1000×2.048=2048 個(gè)數(shù)據(jù)點(diǎn)。 為一隨時(shí)間變化的At隨機(jī)過(guò)程，由于隨機(jī)過(guò)程有很多類型，如維納過(guò)程、正態(tài)隨機(jī)過(guò)程，本文采用了均值為0，方差為 1 的正態(tài)隨機(jī)過(guò)程來(lái)作為演示，來(lái)代替 ，高斯白噪聲采用強(qiáng)度為 2 的高斯At白噪聲代替 ，其帶通濾波后為 。其中濾波器采用的是契比雪夫數(shù)字濾波器。nt dnt可得到 x(t)如下圖所示：0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2-6-4-20246原始輸入信號(hào)1、周期圖法matlab 中的周期圖功率譜法原理是通過(guò)計(jì)算采樣信號(hào)的 FFT，獲得離散點(diǎn)的幅度，再根據(jù)幅度與功率之間的關(guān)系，轉(zhuǎn)換為離散點(diǎn)的功率，再通過(guò)坐標(biāo)變換將離散點(diǎn)的功率圖轉(zhuǎn)換為連續(xù)功率譜密度。Step1:計(jì)算采樣信號(hào) x(n)的 DFT，使用 FFT 方法來(lái)計(jì)算。如果此處將復(fù)頻率處的幅度對(duì)稱到物理實(shí)際頻率，得到的就是單邊譜，否則就是雙邊譜Step2:根據(jù)正余弦信號(hào)功率與幅度的關(guān)系以及直流功率與幅度的關(guān)系，將幅度轉(zhuǎn)換為離散功率譜。Step3:對(duì)橫縱坐標(biāo)進(jìn)行轉(zhuǎn)換，橫坐標(biāo)乘以頻率分辨率轉(zhuǎn)換為實(shí)際連續(xù)物理頻率，縱坐標(biāo)除以頻率分辨率轉(zhuǎn)換為功率譜密度。調(diào)用 MATLAB 中自帶的 matlab 中Pxx,f=periodogram(x,window,nfft,fs)函數(shù)可得計(jì)算結(jié)果如下：0 50 100 150 200 250 300 350 400 450 500f/Hz00.511.52功率/db周期圖法求功率譜2、ARMA 方法參數(shù)模型估計(jì)的思想是:假定研究的過(guò)程 X(n)是一個(gè)輸入序列 u(n)激勵(lì)一個(gè)線性系統(tǒng) H(z)的輸出。有已知的 X(n)，或其自相關(guān)函數(shù)來(lái)估計(jì) H(z)的參數(shù)。由 H(z)的參數(shù)來(lái)估計(jì) X(n)的功率譜。不論 X(n)是確定性信號(hào)還是隨機(jī)信號(hào)，u(n)與 X(n)之間總有如下輸入輸出關(guān)系：10()()()pqkkxnaxnbun0kh對(duì)以上兩個(gè)式子兩邊分別取 Z 變換，并假定 b0=1，可得()BzHA其中 , , 。1()pkkAzaz1()qkkBz0()()kkh為了保證 H(z)是穩(wěn)定的最小相位系統(tǒng)，A(z) 和 B(z)的零點(diǎn)都應(yīng)該在單位圓內(nèi)。假定 u(n)是一個(gè)方差為 的白噪聲序列，由隨機(jī)信號(hào)通過(guò)線性系統(tǒng)的理論可知，輸出序列 X(n)的功2率譜為： 222*()()() jwjwjjwx jBeBePeAAARMA 階數(shù)確定：本題目采用 AIC 準(zhǔn)則確定 ARMA 的階數(shù)。分別計(jì)算 p、q 從 1 到 20 階數(shù)的計(jì)算出AIC（p,q）,如下圖所示，當(dāng)橫坐標(biāo)大概為 230 左右時(shí)，AIC(p,q)取得最小，將此時(shí)的 p,q作為帶入到模型即可。0 50 100 150 200 250 300 350 4000.70.80.911.11.21.3AIC(p,q)ARMA 法譜估計(jì)結(jié)果：0 50 100 150 200 250 300 350 400 450 500f/Hz-5051015202530振幅/dBARMA法 (AIC準(zhǔn)則 )3、Burg 最大熵法Burg 算法的具體實(shí)現(xiàn)步驟：步驟 1 計(jì)算預(yù)測(cè)誤差功率的初始值和前、后向預(yù)測(cè)誤差的初始值，并令 m = 1。210)(NnxP)(gf步驟 2 求反射系數(shù)NmnmngfK1212)()(步驟 3 計(jì)算前向預(yù)測(cè)濾波器系數(shù)),()(11iaiiam,.m步驟 4 計(jì)算預(yù)測(cè)誤差功率 12)(PK步驟 5 計(jì)算濾波器輸出 )()(11ngfnfmmg步驟 6 令 mm+1，并重復(fù)步驟 2 至步驟 5，直到預(yù)測(cè)誤差功率 Pm 不再明顯減小。最后，再利用 Levinson 遞推關(guān)系式估計(jì) AR 參數(shù)，繼而得到功率譜估計(jì)。Burg 最大熵法譜估計(jì)結(jié)果如下圖：0 50 100 150 200 250 300 350 400 450 500f/fs00.511.5功率譜/dBBurg法譜估計(jì)5.附件中表 sheet1 為某地 2008 年 4 月 28 日凌晨 12 點(diǎn)至 2008 年 5 月 4 日凌晨 12 點(diǎn)的電力系統(tǒng)負(fù)荷數(shù)據(jù)，采樣時(shí)間間隔為 1 小時(shí)，利用 Kalman 方法預(yù)測(cè)該地 5 月 5 日的電力系統(tǒng)負(fù)荷，并給出預(yù)測(cè)誤差（5 月 5 日的實(shí)際負(fù)荷數(shù)據(jù)如表 sheet2） 。解：卡爾曼濾波是以最小均方誤差作為估計(jì)的最佳準(zhǔn)則，來(lái)尋求一套遞推估計(jì)的算法，其基本思想是：采用信號(hào)與噪聲的狀態(tài)空間模型，利用前一時(shí)刻地估計(jì)值和現(xiàn)在時(shí)刻的觀測(cè)值來(lái)更新對(duì)狀態(tài)變量的估計(jì)，求得出現(xiàn)時(shí)刻的估計(jì)值。它適合于實(shí)時(shí)處理和計(jì)算機(jī)運(yùn)算?，F(xiàn)設(shè)線性時(shí)變系統(tǒng)的離散狀態(tài)防城和觀測(cè)方程為：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 時(shí)刻的狀態(tài)矢量和觀測(cè)矢量，F(xiàn)(k,k-1) 為狀態(tài)轉(zhuǎn)移矩陣，U(k) 為k 時(shí)刻動(dòng)態(tài)噪聲，T(k,k-1)為系統(tǒng)控制矩陣， H(k)為 k 時(shí)刻觀測(cè)矩陣， N(k)為 k 時(shí)刻觀測(cè)噪聲?？柭鼮V波的算法流程為：1、預(yù)估計(jì)=F(k,k-1)·X(k-1)X(k)2、計(jì)算預(yù)估計(jì)協(xié)方差矩陣=F(k,k-1)×C(k)×F(k,k-1)'+T(k,k-1)×Q(k)×T(k,k-1)'C(k)Q(k)=U(k)×U(k)'3、計(jì)算卡爾曼增益矩陣K(k)= ×H(k)'×H(k)× ×H(k)'+R(k)-1()C(k)R(k)=N(k)×N(k)'4、更新估計(jì)= +K(k)×Y(k)-H(k)× X(k) X()5、計(jì)算更新后估計(jì)協(xié)方差矩陣= I-K(k)×H(k)× ×I-K(k)×H(k)'+K(k)×R(k)×K(k)'C()C()X(k+1) = (k)C(k+1) =6、重復(fù)以上步驟最終可以獲得如下結(jié)果：20 40 60 80 100 120 140 160 180時(shí)間點(diǎn)數(shù)500100015002000250030003500電力系統(tǒng)負(fù)荷使用 Kalman對(duì)電力系統(tǒng)負(fù)荷數(shù)據(jù)進(jìn)行預(yù)測(cè)負(fù)荷真實(shí)值Kalman預(yù)測(cè)值168170172174176178180182184186188190192時(shí)間點(diǎn)數(shù)250030003500電力系統(tǒng)負(fù)荷使用 Kalman對(duì)電力系統(tǒng)負(fù)荷數(shù)據(jù)進(jìn)行預(yù)測(cè)負(fù)荷真實(shí)值Kalman預(yù)測(cè)值168170172174176178180182184186188190192時(shí)間點(diǎn)數(shù)-10001002003005月5日預(yù)測(cè)值與真實(shí)值誤差預(yù)測(cè)值與真實(shí)值之誤差本題將表中的作為觀測(cè)數(shù)據(jù)，圖中橫坐標(biāo)為 1 表示 2008.4.28 日 1 時(shí)刻數(shù)據(jù)，2 表示2008.4.28 日 2 時(shí)刻的數(shù)據(jù)，一次類推，168 表示 2008.5.5 日 1 時(shí)刻的數(shù)據(jù)。從表中可以看出預(yù)測(cè)誤差的最大值為 300。預(yù)測(cè)誤差的大小與代碼中的 R、Q 值得設(shè)置有關(guān)。Q 越大預(yù)測(cè)誤差越小，但是同時(shí)也表明系統(tǒng)內(nèi)的噪聲很大。本題中取得 Q、R 值均為高斯分布的協(xié)方差。代碼見(jiàn)附錄。6.設(shè)某變壓器內(nèi)部短路后，故障電流信號(hào)分解得到下式：y(t) =20e-t+20sin(t+60°)+12sin(2t+45°)+10sin(3t+30°)+6sin(4t+22.5°)+5sin(5t+36°)式中， ， ， 分別利用小波變換、短時(shí)傅里葉變換和維格納威利分=2f0 =30msf0=50Hz布分析故障電流信號(hào)的時(shí)頻特性。解：（1）小波變換：連續(xù)小波變換的定義： *, 1(,)()()()f us tuCWTusftftds計(jì)算連續(xù)時(shí)間小波變換的 4 個(gè)步驟：1、選取一個(gè)小波，然后將其和待分析信號(hào)從起點(diǎn)開(kāi)始的一部分進(jìn)行相乘積分。 2、計(jì)算相關(guān)系數(shù) c。3、將小波向右移，重復(fù) 1 和 2 的步驟直到分析完整個(gè)信號(hào)。4、將小波進(jìn)行尺度伸縮后再重復(fù) 1，2，3 步驟，直至完成所有尺度的分析。（2）短時(shí)傅里葉變換短時(shí)傅里葉變換定義如下： ,(,)()()itf uSTFuftgftgued()1, 2it iuf tedf （3）維格納威利分布變換維格納威利分布定義如下： de),(WDj-*txtx在 MATLAB 中沒(méi)有維格納威利分布變換的相關(guān)函數(shù)，需要安裝一個(gè) MATLAB 版本的時(shí)頻分析工具箱。調(diào)用里面的函數(shù)即可。小波變換和短時(shí)傅里葉變換 MATLAB 均自帶了相關(guān)的函數(shù)。程序見(jiàn)附錄。代碼運(yùn)行結(jié)果結(jié)果如下：0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1時(shí)間 t/s-20020406080幅值原始信號(hào)小波時(shí)頻圖0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1時(shí)間 t/s050100150200250300頻率f/Hz102030405060700.2 0.3 0.4 0.5 0.6 0.7 0.8時(shí)間 /s050100150200250300350400450500頻率/Hz短時(shí)傅里葉變換結(jié)果Wigner-Ville time-frequency distribution0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1時(shí)間 t/s050100150200250300頻率f/Hz7.假定一電力系統(tǒng)諧波與間諧波信號(hào)的函數(shù)表達(dá)式如下： 12 3 4.1cos(2)cos(50).1cos(50).2cos50ynnnnn其中，采樣頻率為 ，相位 為獨(dú)立的均勻分布 ； 為一噪聲信號(hào)，04Hz14:,U信噪比取為 。分別采用三種現(xiàn)代信號(hào)處理方法進(jìn)行諧波與間諧波頻率提取與譜估計(jì)。dB解：本題目采用的頻率提取的三種方法為小波變換、短時(shí)傅里葉變換和維格納威利分布。采用周期圖法、MUSIC 法、Burg 法進(jìn)行譜估計(jì)。確定出諧波的頻率為 50Hz 和150Hz。0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2時(shí)間 (t/秒 )-3-2-10123幅值原始信號(hào)小波時(shí)頻圖0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2時(shí)間 t/s100200300400500頻率f/Hz0.511.522.533.50.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8時(shí)間 /s0100200300400500頻率/Hz短時(shí)傅里葉變換結(jié)果Wigner-Ville time-frequency distribution0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2時(shí)間 t/s050100150200250300350400450500頻率f/Hz0 50 100150200250300350400450500550600f/Hz00.10.20.30.40.5功率/db周期圖法求功率譜0 50 100 150 200 250 300 350 400 450 500頻率 (f/Hz)-40-200204060功率(dB)MUSIC方法0 50 100150200250300350400450500550600f/fs00.0050.010.015功率譜/dBBurg法譜估計(jì)附錄代碼：第四題：clc;clear;fs=1000;%采樣頻率T=2.048;%采樣時(shí)間t=0:1/fs:T;A = normrnd(0,1,1,length(t);%方差為 1，均值為 0 的高斯分布N=wgn(1,length(t),2);%強(qiáng)度為 2 的高斯白噪聲Dn=bandp(N,390,410,200,450,0.1,30,fs);figure(1);subplot(211);plot(t,N);title('原始高斯白噪聲');subplot(212);plot(t,Dn);title('帶通濾波后高斯白噪聲');Sig=sin(2*pi*100.*t)+1.5*sin(2*pi*300.*t)+A.*sin(2*pi*200.*t)+Dn+N;figure(2);plot(t,Sig);title('原始輸入信號(hào)');axis(0 2.1 -7 7);% 周期圖譜Pxx,f=periodogram(Sig,length(t),fs);%周期圖法figure(3);plot(f,Pxx);title('周期圖法求功率譜');xlabel('f/Hz'); ylabel('功率/db');% ARMA 譜估計(jì) z=iddata(Sig');%將信號(hào)轉(zhuǎn)化為 matlab 接受的格式 RecordAIC=; for p=1:20 %自回歸對(duì)應(yīng) PACF，給定滯后長(zhǎng)度上限 p 和 q for q=1:20%移動(dòng)平均對(duì)應(yīng) ACF m=armax(z(1:length(t),p,q); AIC = aic(m); %armax(p,q)選擇對(duì)應(yīng) FPE 最小，AIC 值最小模型 RecordAIC=RecordAIC;p q AIC; end end for k=1:size(RecordAIC,1) if RecordAIC(k,3)=min(RecordAIC(:,3) %選擇 AIC 最小模型 pa_AIC=RecordAIC(k,1); qa_AIC=RecordAIC(k,2); break; endend mAIC=armax(z(1:length(t),pa_AIC,qa_AIC);Pxx2,f2=freqz(mAIC.c,mAIC.a,fs); P2=(abs(Pxx2).*1).2; P2tol=10*log10(P2); figure(4); plot(f2/pi*fs/2,P2tol); title('ARMA 法(AIC 準(zhǔn)則)');xlabel('f/Hz');ylabel('振幅 /dB'); plot(RecordAIC(:,3);ylabel('AIC(p,q)');% burg 法計(jì)算Pxx,F = pburg(Sig,60,length(t),fs);%burg 法figure(6);plot(F,Pxx);title('Burg 法譜估計(jì)');xlabel('f/fs'); %X 軸坐標(biāo)名稱ylabel('功率譜/dB'); %Y 軸坐標(biāo)名稱% function y=bandp(x,f1,f3,fsl,fsh,rp,rs,Fs)%帶通濾波%使用注意事項(xiàng)：通帶或阻帶的截止頻率與采樣率的選取范圍是不能超過(guò)采樣率的一半%即，f1,f3,fs1,fsh,的值小于 Fs/2%x:需要帶通濾波的序列% f 1：通帶左邊界% f 3：通帶右邊界% fs1：衰減截止左邊界% fsh：衰變截止右邊界%rp：邊帶區(qū)衰減 DB 數(shù)設(shè)置%rs：截止區(qū)衰減 DB 數(shù)設(shè)置%FS：序列 x 的采樣頻率% f1=300;f3=500;%通帶截止頻率上下限% fsl=200;fsh=600;%阻帶截止頻率上下限% rp=0.1;rs=30;%通帶邊衰減 DB 值和阻帶邊衰減 DB 值% Fs=2000;%采樣率%wp1=2*pi*f1/Fs;wp3=2*pi*f3/Fs;wsl=2*pi*fsl/Fs;wsh=2*pi*fsh/Fs;wp=wp1 wp3;ws=wsl wsh;% 設(shè)計(jì)切比雪夫?yàn)V波器；n,wn=cheb1ord(ws/pi,wp/pi,rp,rs);bz1,az1=cheby1(n,rp,wp/pi);%查看設(shè)計(jì)濾波器的曲線h,w=freqz(bz1,az1,256,Fs);h=20*log10(abs(h);y=filter(bz1,az1,x);end第 5 題%本題目需要提醒一點(diǎn)：給的數(shù)據(jù)為觀測(cè)數(shù)據(jù) Z 而不是 Xclc;clear;x1=xlsread('./負(fù)荷數(shù)據(jù) .xls','sheet1');x1=x1(:,2);x2=xlsread('./負(fù)荷數(shù)據(jù) .xls','sheet2');x2=x2(:,2);x=x1;x2;N1=length(x1);N=length(x);A=1;B=0;H=1;w=normrnd(0,1000,1,N);%這里隨便取值v=normrnd(0,1000,1,N); P(1)=16;%隨便取值Z=x;X(1)=24;%隨便取值R=cov(v);Q=cov(w);for i=2:Ntempx=A*X(i-1);%+B*u(i);TempP=A*P(i-1)*A'+Q;K(i)=TempP*H'*1/(H*TempP*H'+R);X(i)=X(i-1)+K(i)*(Z(i)-tempx);P(i)=(1-K(i)*H)*TempP;endt=1:length(Z);figure;plot(t,Z,'b',t,X(t),'r');title('使用 Kalman 對(duì)電力系統(tǒng)負(fù)荷數(shù)據(jù)進(jìn)行預(yù)測(cè)');xlabel('時(shí)間點(diǎn)數(shù)');ylabel('電力系統(tǒng)負(fù)荷');axis tight;legend('負(fù)荷真實(shí)值 ','Kalman 預(yù)測(cè)值');figure;subplot(2,1,1);t=length(x1):length(x);plot(t,x(t),'b',t,X(t),'r');title('使用 Kalman 對(duì)電力系統(tǒng)負(fù)荷數(shù)據(jù)進(jìn)行預(yù)測(cè)');xlabel('時(shí)間點(diǎn)數(shù)');ylabel('電力系統(tǒng)負(fù)荷');axis tight;legend('負(fù)荷真實(shí)值','Kalman 預(yù)測(cè)值');set(gca,'XTick',length(x1):2:length(x);subplot(2,1,2);error=Z-X'plot(t,error(t);title('預(yù)測(cè)值與真實(shí)值之誤差');xlabel('時(shí)間點(diǎn)數(shù)');set(gca,'XTick',length(x1):2:length(x);ylabel('5 月 5 日預(yù)測(cè)值與真實(shí)值誤差');axis tight;第六題：% 小波變換 clc;clear; close all;f=50;%信號(hào)頻率 oumiga=2*pi*f;N_sample=2048;%總采樣點(diǎn)數(shù) Fs=1000;%采樣頻率 t=0:1/Fs:1;Tao=0.03;A=1;%信號(hào)幅度 x = 20*exp(-t/Tao)+20*sin(oumiga*t+pi/3)+12*sin(2*oumiga*t+pi/4)+10*sin(3*oumiga*t+pi/6)+6*sin(4*oumiga*t+pi/8)+5*sin(5*oumiga*t+pi/5); % 信號(hào)函數(shù)表達(dá)式figure;plot(t,x);title('原始信號(hào) ');xlabel('時(shí)間 t/s','FontSize',14);ylabel('幅值','FontSize',14);%原信號(hào)函數(shù)wavename='cmor3-3'totalscal=256;Fc=centfrq(wavename); %小波中心頻率c=2*Fc*totalscal;scals=c./(1:totalscal);f=scal2frq(scals,wavename,1/Fs); % 將尺度轉(zhuǎn)換為頻率coefs=cwt(x,scals,wavename); % 求連續(xù)小波系數(shù)figure;imagesc(t,f,abs(coefs); colorbar;xlabel('時(shí)間 t/s','FontSize',14);ylabel('頻率 f/Hz','FontSize',14);title('小波時(shí)頻圖','FontSize',16);axis(0 1 0 300);% 短時(shí)傅里葉變換S,F,T,P=spectrogram(x,256,250,256,Fs);figure;surf(T,F,10*log10(P),'edgecolor','none'); axis tight;view(0,90);xlabel('時(shí)間/s'); ylabel(' 頻率/Hz');title('短時(shí)傅里葉變換結(jié)果');% Wigner-Ville time-frequency distribution.X=hilbert(x');tfr,t,f=tfrwv(X);figure;contour(t/Fs,f*Fs,abs(tfr);xlabel('時(shí)間 t/s');ylabel('頻率 f/Hz');title('Wigner-Ville time-frequency distribution');axis(0 1 0 300)%第七題：clc;clear;close all;% 參數(shù)設(shè)置Fs = 1024; %采樣頻率n = 0:1/Fs:2.01;%采樣時(shí)間N = length(n); % 采樣點(diǎn)W1=0.001*cos(2*pi*n*10+unifrnd(-pi,pi)+cos(2*pi*50*n+unifrnd(-pi,pi)+0.1*cos(2*pi*n*150+unifrnd(-pi,pi)+0.002*cos(2*pi*n*50+unifrnd(-pi,pi);% 原始信號(hào)x1=awgn(W1,20); %加入噪聲%原信號(hào)輸出figure;plot(n,x1);xlabel('時(shí)間(t/秒)','FontSize',10);ylabel('幅值','FontSize',10); axis(0 2.05 -3 3);title('原始信號(hào) ');% 小波變換wavename='cmor3-3'totalscal=256;Fc=centfrq(wavename); %小波中心頻率c=2*Fc*totalscal;scals=c./(1:totalscal);f=scal2frq(scals,wavename,1/Fs); % 將尺度轉(zhuǎn)換為頻率coefs=cwt(x1,scals,wavename); % 求連續(xù)小波系數(shù)figure;imagesc(n,f,abs(coefs); colorbar;xlabel('時(shí)間 t/s','FontSize',14);ylabel('頻率 f/Hz','FontSize',14);title('小波時(shí)頻圖','FontSize',16);% 短時(shí)傅里葉變換S,F,T,P=spectrogram(x1,256,250,256,Fs);figure;surf(T,F,10*log10(P),'edgecolor','none'); axis tight;view(0,90);xlabel('時(shí)間/s'); ylabel(' 頻率/Hz');title('短時(shí)傅里葉變換結(jié)果');% 維格納威利分布X=hilbert(x1');tfr,t,f=tfrwv(X);figure;contour(t/Fs,f*Fs,abs(tfr);xlabel('時(shí)間 t/s');ylabel('頻率 f/Hz');title('Wigner-Villetime-frequency distribution');% 周期圖譜估計(jì)Pxx,f=periodogram(x1,length(x1),Fs);%周期圖法figure;plot(f,Pxx);title('周期圖法求功率譜');xlabel('f/Hz'); ylabel('功率/db');set(gca,'XTick',0:50:600);% MUSIC 方法 譜估計(jì)nfft=1024;figure;pmusic(x1,7,1.1,nfft,Fs,32,16); grid on;xlabel('頻率(f/Hz)','FontSize',10);ylabel('功率(dB)','FontSize',10);title('MUSIC 方法');% burg 法 譜估計(jì)Pxx,F = pburg(x1,60,length(x1),Fs);%burg法figure;plot(F,Pxx);title('Burg 法譜估計(jì)');xlabel('f/fs'); %X 軸坐標(biāo)名稱ylabel('功率譜/dB'); %Y 軸坐標(biāo)名稱set(gca,'XTick',0:50:600);