程序.txt

这个是短时傅立叶变换程序,希望对大家有用了!

%高通滤波
fs=2000
M=256;
n=[0:1:M-1];
x1=sin(2*pi/fs*100*n);
x2=sin(2*pi/fs*800*n);
x=x1+x2;
plot(x)
plot(abs(fft(x)))
Fs=1000;
wp=300*2/Fs;
ws=200*2/Fs;
Rp=1;
Rs=15;
Nn=128
[N,Wn]=buttord(wp,ws,Rp,Rs)
[b,a]=butter(N,Wn,'high');
freqz(b,a,Nn,fs);
y=filter(b,a,x);
plot(abs(fft(y)))
plot(y)



%fft
fs=1000;
N=1024;
n=0:N-1;
t=n/fs;
f1=50;f2=120;
x=sin(2*pi*f1*t)+sin(2*pi*f2*t);
x=x+2*randn(1,length(t));
y=fft(x,N);
mag=abs(y);
f=(0:length(y)-1)'*fs/length(y);
subplot(221)
plot(f(1:N/2),mag(1:N/2));
xlabel('Frequency(Hz)');
ylabel('Magnitude');
title('N=1024 With Noise')



 


 z=amgauss(128,50,40).*fmhyp(128,[1 0.5],[32 0.1]);


%STFT
k=4;T=5;
fc=k*T;         %20
fs=4*fc;        %80
Ts=1/fs;        
N=T/Ts;         %400
x=zeros(1,N);
t=0:N-1;
x=exp(j*k*pi*(t*Ts).^2);
subplot(221);
plot(t*Ts,real(x));
X=fft(x);
X=fftshift(X);
subplot(222);
plot((t-N/2)*fs/N,abs(X));
Nw=20;
L=Nw/2;
Tn=(N-Nw)/L+1;      %39
nfft=32;
TF=zeros(Tn,nfft);
for i=1:Tn
b=hamming(2L);
xw=x((i-1)*L+1:i*L+L).*b';
temp=fft(xw,nfft);
temp=fftshift(temp);
TF(i,:)=temp;
end
subplot(223);
fnew=((1:nfft)-nfft/2)*fs/nfft;
tnew=(1:Tn)*L*Ts;
[F,T]=meshgrid(fnew,tnew);
mesh(F,T,abs(TF));
subplot(224);
contour(F,T,abs(TF));



k=4;T=5;
fc=k*T;         %20
fs=4*fc;        %80
Ts=1/fs;        
N=T/Ts;         %400
x=zeros(1,N);
t=0:N-1;
x=amgauss(400,50,40).*fmlin(400,0.05,0.5,50);
x=x'
subplot(221);
plot(t*Ts,real(x));
X=fft(x);
X=fftshift(X);
subplot(222);
plot((t-N/2)*fs/N,abs(X));
Nw=20;
L=Nw/2;
Tn=(N-Nw)/L+1;      %39
nfft=32;
TF=zeros(Tn,nfft);
for i=1:Tn
b=hamming(20);
xw=x((i-1)*10+1:i*10+10).*b';
temp=fft(xw,nfft);
temp=fftshift(temp);
TF(i,:)=temp;
end
subplot(223);
fnew=((1:nfft)-nfft/2)*fs/nfft;
tnew=(1:Tn)*L*Ts;
[F,T]=meshgrid(fnew,tnew);
mesh(F,T,abs(TF));
subplot(224);
contour(F,T,abs(TF));



fs=2000;
M=256;
n=[0:1:M-1];
x1=sin(2*pi/fs*100*n);
x2=sin(2*pi/fs*800*n);
x=x1+x2;
subplot(211);
plot(x);
subplot(212);plot(abs(fft(x)));




N=400;t=0:N-1;
x=zeros(1,N);
x(100:200)=sin(2*pi*(t(100:200)-100)/20);
plot(x);
subplot(221);
plot(x);
X=fft(x);
X=fftshift(X);
subplot(222);
plot(abs(X));
Nw=10;
L=Nw/2;
Tn=(N-Nw)/L+1;      %39
nfft=32;
TF=zeros(Tn,nfft);
for i=1:Tn
b=hamming(Nw);xw=x((i-1)*L+1:i*L+L).*b';
temp=fft(xw,nfft);
temp=fftshift(temp);
TF(i,:)=temp;
end
subplot(223);
mesh(abs(TF));
subplot(224);
contour(abs(TF))

N=400;t=0:N-1;
x=zeros(1,N);
x(50:100)=cos(pi*(t(50:100)-50)/10);
x(300:360)=cos(pi*(t(300:360)-300)/5);


fs=800;
t=[0:799];
x1=sin(2*pi/fs*200*t);
x2=sin(2*pi/fs*50*t);
x=x1+x2;
subplot(211);
plot(x);
subplot(212);
plot(abs(fft(x)));






wp=150*2/fs;
ws=100*2/fs;
Rp=1;
Rs=15;
Nn=128
[N,Wn]=buttord(wp,ws,Rp,Rs)
[b,a]=butter(N,Wn,'high');
freqz(b,a,Nn,fs);


y=filter(b,a,x);
subplot(211);
plot(y)
subplot(212);
plot(abs(fft(y)))


N=800;
x=zeros(1,N);
x(50:200)=y(50:200);
subplot(221);
plot(x);
X=fft(x);
X=fftshift(X);
subplot(222);
plot(abs(X));
Nw=40;
L=Nw/2;
Tn=(N-Nw)/L+1;      %39
nfft=32;
TF=zeros(Tn,nfft);
for i=1:Tn
b=hamming(Nw);xw=x((i-1)*L+1:i*L+L).*b';
temp=fft(xw,nfft);
temp=fftshift(temp);
TF(i,:)=temp;
end
subplot(223);
mesh(abs(TF));
subplot(224);
contour(abs(TF))