sanbeipin002.m

此例程用于求任意信号的1/3倍程频谱,是我毕业设计中的一个小小的部分,经验证算法正确.

load clear1.mat cl0;
x=cl0;
sf=5000;
f=[1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.30 8.00];
fc=[f,10*f,100*f,1000*f,10000*f];
oc6=2^(1/6);
nc=length(fc);
n=length(x);
nfft=2^nextpow2(n);
a=fft(x,nfft);
for j=1:nc
    fl=fc(j)/oc6;
    fu=fc(j)*oc6;
    nl=round(fl*nfft/sf+1);
    nu=round(fu*nfft/sf+1);
    if fu>sf/2
       m=j-1;break
    end
    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(j)=sqrt(var(real(b(1:n))));
end
figure (5);
subplot(2,1,1);
t=0:1/sf:(n-1)/sf;
plot(t,x);
title('原始加速度信号','color','b','FontSize',14);
xlabel('时间(s)');
ylabel('加速度(g)');
grid on;
subplot(2,1,2);
plot(fc(1:m),yc(1:m),'r');
title('三分之一倍频程谱','color','b','FontSize',14);
xlabel('频率(HZ)');
ylabel('有效值');
grid on;