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

《西南交通大學信號處理期末作業(yè)》由會員分享，可在線閱讀，更多相關(guān)《西南交通大學信號處理期末作業(yè)（14頁珍藏版）》請在裝配圖網(wǎng)上搜索。

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