那位给个matlab的短时傅立叶变换的程序!

wrx52574800

如题

eight

原帖由 wrx52574800 于 2007-5-20 16:45 发表

                               
登录/注册后可看大图

如题

请移步信号处理方法版块，或直接搜索一下论坛

tujian510

这两天搞毕业设计，编的用Matlab中的fft分析含噪声的信号，其中有直接FFT，时域同步平均后FFT，频域同步平均后FFT。可以参考。

%%This procedure is to investigate the FFT of sine signal, sine signal with
%%noise.  The sine signal with noise is treated by direct FFT,
%%FFT after piece-wise averaging in time domain, and FFT after
%%piece-wise averaging  frequency domain
%%
clear all
pack


w=10;
t=0:(2*pi/w/500):(300*2*pi/w);
x=sin(w*t);

n=length(t);

ss=888;
randn('seed',ss);
z=randn(1,n);
z1=z-mean(z);
z2=z1/max(abs(z1));
z2=z2*500;

figure(1)
plot(t,x)  % The sine signal without noise


fre=(1/(2*pi/w/500))*(1:n)/n;

nn=fix(n/2);
y1=fft(x);  % The FFT of sine signal without noise
y1=(abs(y1)).^2/n;
y2=y1(2:(nn+1));
fre1=2*pi*fre(1:nn);

figure(2)
plot(fre1,y2)

x=x+z2;

y1=fft(x);  % The FFT of sine signal with noise
y1=(abs(y1)).^2/n;
y2=y1(2:(nn+1));

figure(3)
plot(fre1,y2)         %The sine signal with noise

%% Averaging in time domain

M=10000;
xx=zeros(1,M);
m1=1;m2=M;

nnnn=fix(300*500/M);

for nnn=1:nnnn
    xxx=x(m1:m2);
    xx=xx+xxx;
    m1=m1+M;
    m2=m2+M;
end

fre=(1/(2*pi/w/500))*(1:M)/M;

nn=fix(M/2);
xx=xx/nnnn;
y1=fft(xx);
y1=(abs(y1)).^2/M;
y2=y1(2:(nn+1));
fre1=2*pi*fre(1:nn);

figure(4)
plot(fre1,y2)



% Next is averaging in frequency domain
% In fact, the averaging in time domain is identical with that in frequency
% domain

MM=10000;
xx=zeros(1,MM);
m1=1;m2=MM;

nnnn=fix(300*500/MM);

for nnn=1:nnnn
    xxx=x(m1:m2);
    xxx=fft(xxx);
    %xxx=sqrt((abs(xxx)).^2);
    xx=xx+xxx;
    m1=m1+MM;
    m2=m2+MM;
end

y1=xx/nnnn;
y1=(abs(y1)).^2/MM;


fre=(1/(2*pi/w/500))*(1:MM)/MM;

nn=fix(MM/2);
y2=y1(2:(nn+1));
fre1=2*pi*fre(1:nn);

figure(5)
plot(fre1,y2)