1/2倍频程 如何用Matlab求解

logcos

1/2倍频程 如何用Matlab求解                                                                                                                                                                                                                                                                                                                        老师给的要求：
读出一个Wav噪声 得到其原波形图
对声音信号进行傅里叶变换 得到其频谱图
再运行傅里叶变换 得到倒谱
得到声音信号的1/2倍频程图

我的写的程序如下： 不会救1/2倍频程

[y,fs,bits]=wavread('d:\wav.wav');
t=length(y);
subplot(3,1,1)
plot(y);
legend('波形图');
xlabel('时间(s)');
ylabel('幅度');
%---------------
x=fft(y.*hamming(t));
fm=5000*length(x)/fs;
f=(0:fm)*fs/length(x);
subplot(3,1,2);
plot(f,20*log10(abs(x(1:length(f))+eps)));
legend('频谱图');
xlabel('频率Hz');
ylabel('频率幅度db');

c=fft(log(abs(x)+eps));
ms1=fs/1000;
ms20=fs/50;
q=(ms1:ms20)/fs;
subplot(3,1,3);
plot(q,abs(c(ms1:ms20)));
legend('倒谱图');
xlabel('倒频');
ylabel('倒频谱幅度');

请问对不，再请高手指点一下，1/2倍频程怎么求，谢谢～

logcos

[y,Fs,bits]=wavread('d:\wav\wav.wav');
N=3;
Fc=Fs/100;
pi = 3.14159265358979;
f1 = Fc/(2^(1/6));
f2 = Fc*(2^(1/6));
Qr = Fc/(f2-f1);
Qd = (pi/2/N)/(sin(pi/2/N))*Qr;
alpha = (1 + sqrt(1+4*Qd^2))/2/Qd;
W1 = Fc/(Fs/2)/alpha;
W2 = Fc/(Fs/2)*alpha;
plot(butter(N,[W1,W2]));

这样写得出来的图是这样的，不对呀～

w89986581

function [fc,yc] = octspectrum(yy,fs)
%yy信号;fs采样率
f=1/4:1/4:16;
fc=2.^f;
oc6=2^(1/4);
nc=length(fc);
n=length(yy);
nfft=2^nextpow2 ;
a=fft(yy,nfft);
for jj=1:nc,
    fl=fc(jj)/oc6;
    fu=fc(jj)*oc6;
    nl=round(fl*nfft/fs+1);
    nu=round(fu*nfft/fs+1);
    b=zeros(1,nfft);
    b(nl:nu)=a(nl:nu);
    b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
    c=ifft(b,nfft);
    %%%计算对应每个中心频率段的有效值
    yc(jj)=sqrt(var(real(c(1:n))));
end
return