fa1.m

三分之一倍频程表示的短时傅里叶变换

%三分之一倍频程处理
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%打开wav文件，存储单声道声音数据在新的文本文件中
[FID,fs,nbit]=wavread('C:\Documents and Settings\Administrator\桌面\毕设\5s_ideling.wav');
A=FID(:,1);

fid=fopen('lorry_1.txt','w');    %以写的方式打开或建立一个新文件
count=fwrite(fid,A,'double');    %把矩阵A中的数据读入该文件
fid=fopen('lorry_1.txt','r')     %以只读的方式打开文件
x=fread(fid,'double');           %读入时程数据存成列向量
status=fclose(fid);              %关闭数据文件


%定义三分之一倍频程的中心频率
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);
%大于并最接近n的2的幂次方长度
nfft=2^nextpow2(n);
%FFT变换
a=fft(x,nfft);
for j=1:nc
    %下限频率
    fl=fc(j)/oc6;
    %上限频率
    fu=fc(j)*oc6;
    %下限频率对应的序号
    nl=round(fl*nfft/fs+1);
    %上限频率对应的序号
    nu=round(fu*nfft/fs+1);
    %如果上限频率大于折叠频率则循环中断
    if fu>fs/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/fs:(n-1)/fs;
plot(t,x);
xlabel('时间(s)');
ylabel('声压级(dB)');
grid on
%绘制三分之一倍频程有效值图形
subplot(2,1,2);
%plot(fc(1:nc),yc(1:nc));
plot(fc(1:nc),yc(1:nc));
xlabel('频率(Hz)');
ylabel('有效值');
grid on


    
    
%建立离散时间列向量
%t=(0:1/fs:(n-1)/fs)';
%计算趋势项的多项式待定系数向量a
%a=polyfit(t,x,4); % 4为拟合多项式的阶数
%用x减去多项式系数a生成的趋势项