傅里叶变换对噪声进行频谱分析

%%%%%%傅里叶变换/逆变换/短时傅里叶变换%%%%%%

[y,Fs]=wavread('Noise_5l_cai.wav'); %读出信号，采样率和采样位数。
%sound(y,Fs);
y=y(:,1); %我这里假设你的声音是双声道，我只取单声道作分析，如果你想分析另外一个声道，请改成y=y(:,2)
figure(1)
sigLength=length(y); 
%
Nshow=round(sigLength/10); %显示的点数（频率段，因频谱对称只需一半内）
%
t=(0:length(y)-1)/Fs; %计算时间轴
subplot(3,1,1);plot(t,y);xlabel('Time(s)');%在第一个窗口画波形
grid on;

n=0:sigLength-1;
Y= fft(y,length(y)); %在sigLength这个有限区间内做快速傅立叶变换 
%mag=abs(Y);
f=Fs*n/sigLength;
subplot(3,1,2);plot(f(1:Nshow),abs(Y(1:Nshow)));
xlabel('Frequency(Hz)'); 
%在第一个窗口画率谱, 
grid on;

subplot(3,1,3);
xifft=ifft(Y);
%sound(real(xifft),Fs);
plot(t,xifft);xlabel('Time(s)')
grid on;

figure(2)
Nw=200;          %窗函数长 window length
L=Nw/2;         %窗函数每次移动的样点数
Ts=floor((sigLength-Nw)/L)+1;     %计算把数据y共分成多少段 round四舍五入，floor向下取整
nfft=512*5;        % FFT的长度
TF=zeros(Ts,nfft/10);  %将存放三维图谱，先清零 /20-add 只显示有用的频率
for i=1:Ts
    xw=y((i-1)*L+1:i*L+L);   %取一段数据
    temp=fft(xw,nfft);       % FFT变换
    %temp=fftshift(temp);     %频谱以0频为中心
    TF(i,:)=temp(1:nfft/10)';            %把谱图存放在TF中代表TF矩阵中的第i行 /20--add
end
fn=(1:nfft/10)*Fs/nfft; %/20--add 只显示一半的点
tn=(1:Ts)*Nw/2/Fs;
[T,F]=meshgrid(tn,fn);
mesh(T,F,abs(TF.'));  %三维绘图
axis tight;
title('三维时频图'); 
ylabel('Frequency(Hz)'); 
xlabel('Time(s)'); 
grid on;