三分之一倍频程程序(2)

%上限频率对应的序号 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

%绘制输入时程曲线图形 subplot(2,1,1); t=0:1/sf:(n-1)/sf; plot(t,x);

xlabel('时间(s)'); ylabel('加速度(g)'); grid on;

%绘制三分之一倍频程有效值图形 subplot(2,1,2);

plot(fc(1:m),yc(1:m));

xlabel('频率(Hz)'); ylabel('有效值'); grid on;

%保存倍频程数据 fid=fopen(fno,’w’); for k=1:m;

fprintf(fid,’%f %f\\n’,fc(k),yc(k)); end

status=fclose(fid);