Fs bit]=wavread('50.wav'); %讀取波形文件，s為數據，Fs為抽樣頻率，bit為位寬
figure(1) %新建壹個窗口
%小波去噪
[thr,sorh,keepapp]=ddencmp('den','wv',s); %返回小波除噪和壓縮後的信號
[s,cxc,lxc,perf,perfl2]=wdencmp('gbl',s,'db3',3,thr,sorh,keepapp); %通過門限閾值處理得到小波系數
%求頻譜
nfft=length(s); %信號長度，即抽樣點數
p=floor(length(s)/3)+1; %AR預測的階數
[spsd,f]=pyulear(s,p,nfft,Fs,'onesided'); %求功率譜密度
pmax=max(spsd); %功率譜密度最大分量
spsd=spsd/pmax; %歸壹化處理
%顯示
plot(f(1:200),spsd(1:200)); %對應頻率點的功率譜畫圖
a=find(spsd==1) %尋找歸壹化後為1的頻率點的位置
f(a) %輸出最大功率譜密度對應的頻率