这段程序怎么解释？

kingsword1

各位高手，这是本论坛里面的一个帖子里面其中的一段程序，是用来画三维图形的时频分析。程序如下：
%%短时傅立叶变换
%解析信号x=exp(j*pi*k*t.^2)
clear,clc,close all
figure(1)
k=4;T=5;
fc=k*T;
fs=3*fc;
Ts=1/fs;
N=T/Ts;
x=zeros(1,N);
t=0:N-1;
x=exp(j*k*pi*(t*Ts).^2);
% x=awgn(x,-3,'measured');
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;
nfft=32;
TF=zeros(Tn,nfft);
for i=1:Tn
    xw=x((i-1)*10+1:i*10+10);
    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))
其中我对后面画三维图那里的程序不是很明白，比如：NW是什么意思？L又是什么意思，为什么要NW/2？Tn=（N-NW）/L+1又如何解释？nfft的取值有什么说法没？

songzy41

这一段程序中的关键是进行短时傅里叶变换，然而构成一个三维谱图。在取数时每
段数据的窗函数是移动10个数据样点，而窗函数每次取20个数据样点。
Nw=20;              %窗函数长
L=Nw/2;             %窗函数每次移动的样点数
Tn=(N-Nw)/L+1;      %计算把数据x共分成多少段
nfft=32;            %FFT的长度
TF=zeros(Tn,nfft);  %将存放三维谱图，先清零
for i=1:Tn
    xw=x((i-1)*10+1:i*10+10);  %取一段数据
    temp=fft(xw,nfft);         %FFT变换
    temp=fftshift(temp);       %频谱以0频为中心
    TF(i,:)=temp;              %把谱图存放在TF中
end

kingsword1

楼上的高手，先谢谢你给出的解释！现在我要解决的就是我有一段采集的数据，采样频率为500Hz，我想取其中的5秒来分析，不知道我的程序改按照上面的程序怎么改？能帮帮我吗？万分感激！

kingsword1

我试着用上面的程序编了一下，出来个这样的图（附图1），是双边频的，怎么把他改成单边频的呢？还有别人的图是这样的（附图2），这个程序怎么改才能变成图2那样呢？
程序如下：
Nw=20;fs=500;t=5;
Ts=1/fs;N=t/Ts;
L=Nw/2;
Tn=(N-Nw)/L+1;
nfft=1024;
TF=zeros(Tn,nfft);
for i=1:Tn
    xw=x((i-1)*10+1:i*10+10);
    temp=fft(xw,nfft);
    temp=fftshift(temp);
    TF(i,:)=temp;
end
figure(2);
fnew=((1:nfft)-nfft/2)*fs/nfft;
tnew=(1:Tn)*L*Ts;
[F,T]=meshgrid(fnew,tnew);
mesh(T,F,abs(TF));

songzy41

程序作如下的改动，可以只显示单边的：
Nw=20;
L=Nw/2;
Tn=(N-Nw)/L+1;
nfft=32;
TF=zeros(Tn,nfft/2);
for i=1:Tn
    xw=x((i-1)*10+1:i*10+10);
    temp=fft(xw,nfft);
    TF(i,:)=temp(n);
end
fnew=(1:nfft/2)*fs/nfft;
tnew=(1:Tn)*L*Ts;

[F,T]=meshgrid(fnew,tnew);
mesh(F,T,abs(TF))

songzy41

从笫1张图看，有较大的直流分量，在处理前应去除直流分量。

kingsword1

楼上的高手，谢谢你帮我改动程序，但是我运行了你的改动后的程序，出错呀！n没有赋值或者定义，还有就是??? Subscripted assignment dimension mismatch.原来的程序可以出图呀，现在的不行了，能否再修改一下，完善一下你的程序呢，谢谢！

fafan

什么环境运行的？

kingsword1

MATLAB7的环境下呀

kingsword1979

有人可以帮我改改上面的程序吗？

songzy41

原帖由 kingsword1 于 2006-10-29 13:14 发表
楼上的高手，谢谢你帮我改动程序，但是我运行了你的改动后的程序，出错呀！n没有赋值或者定义，还有就是??? Subscripted assignment dimension mismatch.原来的程序可以出图呀，现在的不行了，能否再修改一下，完 ...

把上程序中的语句：
TF(i,:)=temp(n);
改成
TF(i,:)=temp;
便可运行了。

weiyuperfect

非常感谢！我在考虑采用短时傅立叶变换来识别信号，不知道会不会成功。非常感谢

weiyuperfect

1.temp=fftshift(temp);       %频谱以0频为中心
以0频为中心是什么意思呢？4楼的第一个图中为什么频率的值是从-400到400呢？频率值怎么会有负值呢？

2.这里的横坐标是时间，纵坐标是频率，竖坐标是短时傅立叶分析后得到的数据吧？

kexin

回复5楼，
    temp=fft(xw,nfft);应改为：
temp=fft(xw,nfft/2);这样就可以正确运行了

xianiao

plot((t-N/2)*fs/N,abs(X))这句什么意思，为什么这样画图？
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(T,F,abs(TF))这些是什么？对这些图还不懂，希望高人能给解释下！谢谢！