求短时傅立叶变换的程序包

weiyuperfect

请问各位大侠，哪位有短时傅立叶变换的程序？非常感谢！！

songzy41

在时频分析处理工具箱(tftb)中带有短时傅立叶变换的程序。

weiyuperfect

谢谢您！
1.您说的tftb中哪个是短时傅立叶变换的程序呢？函数名是什么？是不是能直接用？如何调用呢？

2.我算出了一三层框架一层加速度信号，即每隔0.02秒对应一个加速度值，想对这个信号进行短时傅立叶变换，不知道该怎么做。是否可以用下面的程序呢？可我觉得下面的程序并没有把每次得到的傅立叶变换的数据都表示出来啊？应该怎么做呢？

3.另外我在论坛里看了下面的程序：
>> fs=3200;     %采样频率
>> t=[0:1/fs:0.1-1/fs];
>> a=1.414*ones(size(t));
>> index=find(0.035<=t&t<=0.075);
>> a(index)=1.414*0.7;
>> y=a.*sin(100*pi*t);
>> subplot(211);
plot(t,y);
xlabel('时间');
ylabel('幅值');
>> title('电压暂降信号 时域波形');
>> %短时傅立叶变换
>> N=length(y);
>> Nw=64;                     %窗函数长
>> L=16;                       %窗函数每次移动的样点数
>> Tn=(N-Nw)/L+1;              %计算把数据x共分成多少段
>> nfft=64;                   %FFT的长度
>> TF=zeros(Tn,nfft);          %将存放三维谱图，先清零
>> for i=1:Tn
n1=L*(i-1)+1;%傅立叶变换起点，
n2=L*(i-1)+Nw;%傅立叶变换终点
sf=fft(y(n1:n2),nfft);     %FFT变换
TF(i,:)=sf;                %把谱图存放在TF中
end
>> subplot(212);
>> fnew=((1:nfft)-nfft/2)*fs/nfft;
>> tnew=(1:Tn)*L/fs;
[F,T]=meshgrid(fnew,tnew);
>> plot(tnew,2*abs(TF(:,2))/(sqrt(2)*Nw),'b');
>> xlabel('时间');
>> ylabel('有效值');
>> title('短时傅立叶跟踪电压暂降信号的有效值波形');
>> grid on;

请问>> TF=zeros(Tn,nfft);          %将存放三维谱图，先清零
这里的三维谱图三个坐标分别是什么呢?傅立叶变换sf=fft(y(n1:n2),nfft);  得到的数据只是一列的数据啊，其余两个坐标分别是什么呢？如果想画出三维谱图应该用什么命令呢？

4.另外这里的窗函数窗口长和宽都是什么呢？如果您能给予指点我将万分感谢！！！

songzy41

1，在tftb中有tfrstft函数，便是做STFT的函数，具体用法用help看一下；
2，3中的程序，只能说是计算STFT的一种方法，用tfrstft函数可能更好；
3，TF是一个3维谱，如果用：
imagesc(tnew,fnew(32:64),abs(TF(:,1:33)));
axis('xy');
便能画出3维谱图了。但原程序是取64点的FFT，点数较少，所以谱图较粗；
4，该程序取了nfft=64，便是窗长。

weiyuperfect

这个信号是G1，是由以下程序得到的:
                                             
M=[2573*10^3 0 0;0 2558*10^3 0;0 0 562*10^3];
>> k1=554*10^3;
>> k2=926*10^3;
>> k3=836*10^3;
>> K=[k1+k2 -k2 0;-k2 k2+k3 -k3;0 -k3 k3];
>> M1=inv(M);
>> Q=K*M1;
>> W=sqrt(eig(Q));
>> w1=W(1,1);w2=W(2,1);
>> xishu=2*0.05/(w1+w2)*[w1*w2;1];
>> a=xishu(1,1);b=xishu(2,1);
>> C=a*M+b*K;
>> U=[0;0;0];U1=[0;0;0];U2=[0;0;0];
t=0;
g=[0;0;0];
dt=0.02;sei=1.4;tt=sei*dt;
i=1;
G1=zeros(800,1); G2=zeros(800,1); G3=zeros(800,1);
G=zeros(801,1);
while t<8
dgtt=(sin(t+tt)-sin(t))*[1;1;1];
dptt=-M*dgtt;
dgdt=(sin(t+dt)-sin(t))*[1;1;1];
dpdt=-M*dgdt;
K1=K+6*M/(tt^2)+3*C/tt;
P1=dptt+M*(6*U1/tt+3*U2)+C*(3*U1+0.5*U2*tt);
dutt=inv(K1)*P1;
du2tt=6*(dutt-U1*tt-0.5*U2*tt^2)/tt^2;
du1tt=3*(dutt-U1*tt-(U2*tt^2)/6)/tt;
du2dt=du2tt/sei;
du1dt=U2*dt+0.5*du2dt*dt;
dudt=U1*dt+0.5*U2*dt^2+(du2dt*dt^2)/6;
U=U+dudt;
U1=U1+du1dt;
g=g+dgdt;
G(i,1)=t;
i=i+1;
t=t+dt;
U2=-M1*(M*g+C*U1+K*U);
G1(i,1)=U2(1,1);
G2(i,1)=U2(2,1);
G3(i,1)=U2(3,1);
end
k1=398*10^3;
K=[k1+k2 -k2 0;-k2 k2+k3 -k3;0 -k3 k3];
Q=K*M1;
W=sqrt(eig(Q));
w1=W(1,1);w2=W(2,1);
xishu=2*0.05/(w1+w2)*[w1*w2;1];
a=xishu(1,1);b=xishu(2,1);
C=a*M+b*K;
while t<=16
dgtt=(sin(t+tt)-sin(t))*[1;1;1];
dptt=-M*dgtt;
dgdt=(sin(t+dt)-sin(t))*[1;1;1];
dpdt=-M*dgdt;
K1=K+6*M/(tt^2)+3*C/tt;
P1=dptt+M*(6*U1/tt+3*U2)+C*(3*U1+0.5*U2*tt);
dutt=inv(K1)*P1;
du2tt=6*(dutt-U1*tt-0.5*U2*tt^2)/tt^2;
du1tt=3*(dutt-U1*tt-(U2*tt^2)/6)/tt;
du2dt=du2tt/sei;
du1dt=U2*dt+0.5*du2dt*dt;
dudt=U1*dt+0.5*U2*dt^2+(du2dt*dt^2)/6;
U=U+dudt;
U1=U1+du1dt;
g=g+dgdt;
G(i,1)=t;
t=t+dt;
G1(i,1)=U2(1,1);
G2(i,1)=U2(2,1);
G3(i,1)=U2(3,1);
i=i+1;
U2=-M1*(M*g+C*U1+K*U);
end
plot(G,G

weiyuperfect

非常感谢songzy41！！您回复的帖子我看了非常感激!!!!!!you are a good man!!!

非常抱歉，上面我发的帖子没有写完全，点错了，直接就发出去了，请大家原谅我的马虎失误，还请大家帮忙看看我的问题，多谢了！！

这个信号是G1，是由以下程序得到的:
                                             
M=[2573*10^3 0 0;0 2558*10^3 0;0 0 562*10^3];
>> k1=554*10^3;
>> k2=926*10^3;
>> k3=836*10^3;
>> K=[k1+k2 -k2 0;-k2 k2+k3 -k3;0 -k3 k3];
>> M1=inv(M);
>> Q=K*M1;
>> W=sqrt(eig(Q));
>> w1=W(1,1);w2=W(2,1);
>> xishu=2*0.05/(w1+w2)*[w1*w2;1];
>> a=xishu(1,1);b=xishu(2,1);
>> C=a*M+b*K;
>> U=[0;0;0];U1=[0;0;0];U2=[0;0;0];
t=0;
g=[0;0;0];
dt=0.02;sei=1.4;tt=sei*dt;
i=1;
G1=zeros(800,1); G2=zeros(800,1); G3=zeros(800,1);
G=zeros(801,1);
while t<8
dgtt=(sin(t+tt)-sin(t))*[1;1;1];
dptt=-M*dgtt;
dgdt=(sin(t+dt)-sin(t))*[1;1;1];
dpdt=-M*dgdt;
K1=K+6*M/(tt^2)+3*C/tt;
P1=dptt+M*(6*U1/tt+3*U2)+C*(3*U1+0.5*U2*tt);
dutt=inv(K1)*P1;
du2tt=6*(dutt-U1*tt-0.5*U2*tt^2)/tt^2;
du1tt=3*(dutt-U1*tt-(U2*tt^2)/6)/tt;
du2dt=du2tt/sei;
du1dt=U2*dt+0.5*du2dt*dt;
dudt=U1*dt+0.5*U2*dt^2+(du2dt*dt^2)/6;
U=U+dudt;
U1=U1+du1dt;
g=g+dgdt;
G(i,1)=t;
i=i+1;
t=t+dt;
U2=-M1*(M*g+C*U1+K*U);
G1(i,1)=U2(1,1);
G2(i,1)=U2(2,1);
G3(i,1)=U2(3,1);
end
k1=398*10^3;
K=[k1+k2 -k2 0;-k2 k2+k3 -k3;0 -k3 k3];
Q=K*M1;
W=sqrt(eig(Q));
w1=W(1,1);w2=W(2,1);
xishu=2*0.05/(w1+w2)*[w1*w2;1];
a=xishu(1,1);b=xishu(2,1);
C=a*M+b*K;
while t<=16
dgtt=(sin(t+tt)-sin(t))*[1;1;1];
dptt=-M*dgtt;
dgdt=(sin(t+dt)-sin(t))*[1;1;1];
dpdt=-M*dgdt;
K1=K+6*M/(tt^2)+3*C/tt;
P1=dptt+M*(6*U1/tt+3*U2)+C*(3*U1+0.5*U2*tt);
dutt=inv(K1)*P1;
du2tt=6*(dutt-U1*tt-0.5*U2*tt^2)/tt^2;
du1tt=3*(dutt-U1*tt-(U2*tt^2)/6)/tt;
du2dt=du2tt/sei;
du1dt=U2*dt+0.5*du2dt*dt;
dudt=U1*dt+0.5*U2*dt^2+(du2dt*dt^2)/6;
U=U+dudt;
U1=U1+du1dt;
g=g+dgdt;
G(i,1)=t;
t=t+dt;
G1(i,1)=U2(1,1);
G2(i,1)=U2(2,1);
G3(i,1)=U2(3,1);
i=i+1;
U2=-M1*(M*g+C*U1+K*U);
end

plot(G,G1)就可以得到这个信号的图形，这里我没法把图给传过来，麻烦大侠们运行一下我的程序，就可以得到G1的图形了，现在我用短时傅立叶变换对这个信号分解，程序如下：
X=G1;
T=1:length(X);
N=length(X);
H=20;
TRACE=0;
[TFR,T,F]=TFRSTFT(X,T,N,H,TRACE);
mesh(T,F,TFR)
但看不到有任何突变，我的目的是证明这个信号在8秒时有明显突变，
是不是短时傅立叶变换做不到这个啊?还请您多费心帮忙看看好么？谢谢！！

songzy41

原帖由 weiyuperfect 于 2008-7-30 21:31 发表
但看不到有任何突变，我的目的是证明这个信号在8秒时有明显突变，
是不是短时傅立叶变换做不到这个啊?还请您多费心帮忙看看好么？谢谢！！

用短时傅里叶变换来做奇异点的检测，不是一个好方法。因为在STFT中主要是做谱分析，谱分析就要取一个时间窗，一帧谱是在该时间窗中的“平均值”，时间窗取得短，时间的跟随性会稍好一些，但可能分析不了信号中的低频信号；时间窗取得长，更不能反映出瞬态的特性。例如在3楼给出的例子，是一个电压暂降的情况，电压暂降发生在0.035秒处，即在0.035秒处发生了突变，但分析出来看到是从0.02至0.04秒之间发生的逐步变化，而看不到突变的情况。

weiyuperfect

非常感谢您的帮助，我学到了很多知识，对您非常仰慕！我还有两个问题：
1.频谱以0频为中心是什么意思呢？为什么频率的值是从-400到400呢？频率值怎么会有负值呢？

2.短时傅立叶变换后得到的三维谱图横坐标是时间，纵坐标是频率，竖坐标是短时傅立叶分析后得到的数据吧？

songzy41

原帖由 weiyuperfect 于 2008-7-31 21:41 发表
频率值怎么会有负值呢？
非常感谢您的帮助，我学到了很多知识，对您非常仰慕！我还有两个问题：
1.频谱以0频为中心是什么意思呢？为什么频率的值是从-400到400呢？频率值怎么会有负值呢？

2.短时傅立叶变换后得到的三维谱图横坐标是时间，纵坐标是频率，竖坐标是短时傅立叶分析后得到的数据吧？

1，对于一个余弦信号：x(t)=cos(wt)=[exp(-jwt)+exp(jwt)]/2
其中就包含了正频率部分和负频率部分，傅里叶变换及DFT都分析了正负频率的部分。一般FFT后正频率部分在前（在下标为0~N/2+1之间）；负频率部分在后（在下标为N/2+2~N之间），但用了fftshift后，把0频率调整到中间，就你楼主讲的从负频率到正频率的排列。

2，对的。

weiyuperfect

如题，谢谢您的帮助,我学到了挺多知识的，并对您表示敬意！！

lotus86

你好，我想问一下，你用的MATLAB,是哪个版本啊？我的是06B版，怎么没有tfrstft，这个函数啊:@o
谢谢阿

ChaChing

这不是matlab原有的函数, 需去下载的! (时频工具包)

467744980

用specgram 分析信号怎么不显示负频部分啊？  该怎样让它显示负频？  请指教