short time fourier transform

pengzk · 发表于 2007-4-19 04:51

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

function [Spec] = STFT(Sig,N,WinLen);

% Calculating the Short Time Fourier Transform
%
%    Sig    : the signal to be analyzed
%       N    : the number of frequency bins (default : length(Sig)).
% WinLen    : the length of window used to locate the signal in time.
%
% June 1, 2003
%

if (nargin < 1),
error('At least one parameter required!');
end;

Sig = real(Sig);
SigLen = length(Sig);

if (nargin < 3),
WinLen = SigLen / 4;
end

if (nargin < 2),
N = SigLen;
end

if(N > 512),
   N = 512;
end

WinLen = ceil(WinLen / 2) * 2;
WinFun = exp(log(0.005) * linspace(-1,1,WinLen)'.^2 );
WinFun = WinFun / norm(WinFun);

Lh = (WinLen - 1)/2;
Spec = zeros(N,SigLen) ;

wait = waitbar(0,'Under calculation, please wait...');
for iLoop = 1:SigLen,

waitbar(iLoop/SigLen,wait);
tau = -min([round(N/2)-1,Lh,iLoop-1]):min([round(N/2)-1,Lh,SigLen-iLoop]);
Index = floor(rem(N + tau,N) + 1);
temp = floor(iLoop + tau);
temp1 = floor(Lh+1+tau);
Spec(Index,iLoop) = Sig(temp) .* conj(WinFun(temp1));

end;

close(wait);

Spec = fft(Spec);
Spec = abs(Spec);
f = [0:N/2-1 -N/2:-1]'/N;
t = 1: SigLen;

fmax = 0.5;
fmin = 0.0;

clf
%=====================================================%
% Plot the result                                  %
%=====================================================%
set(gcf,'Position',[20 100 800 500]);
set(gcf,'Color','w');

axes('Position',[0.07 0.37 0.4 0.55]);
mesh(t,f,Spec);
axis([min(t) max(t) fmin fmax]);
ylabel('Frequency');
title('Short Time Fourier Transform');

%=======================================%
% Plot the signal in the local zone    %
%=======================================%
axes('Position',[0.07 0.08 0.4 0.2]);
plot(t,real(Sig),'r');
xlabel('Time');
ylabel('Amp');
axis([1 SigLen min(real(Sig)) max(real(Sig))]);
title('Orignal signal');
grid

%===================================================%
% plot the FFT of the signal                      %
%===================================================%
sptr = fft(Sig);                                     %sptr is the FFT spectrum of the signal
sptrLen = length(sptr);
sptr = sptr(1:round(sptrLen/2)+1);
f = linspace(0,0.5,round(sptrLen/2)+1);
axes('Position',[0.55 0.37 0.4 0.55]);
plot(f,abs(sptr),'b');
axis([min(f) max(f) 0 max(abs(sptr))]);
title('Discrete Fourier Transform');
xlabel('Frequency');
ylabel('Amp');
grid

pengzk · 发表于 2007-4-19 04:52

SampFreq = 30;
t=0:1/SampFreq:40;
sig = sin(12*pi*t);
sig(1:end/2) = sig(1:end/2) + sin(6*pi*t(1:end/2));
sig(end/2+1:end) = sig(end/2+1:end) + sin(18*pi*t(end/2+1:end));
STFT(sig',512,512);

zhlong · 发表于 2007-4-19 09:35

多谢共享自己的成果

songzy41 · 发表于 2007-4-19 15:48

感谢共享，不过很遗憾的是STFT只能处理实数信号，因为在STFT函数中有
Sig = real(Sig);
语句。如果能处理复数信号则更完善了。

ghb063 · 发表于 2008-4-1 22:36

非常感谢

jerry0204103 · 发表于 2015-1-28 17:20

??? Attempt to execute SCRIPT STFT as a function.

Error in ==> STFTcs at 6
STFT('sig',512,512);

运行出现上边的情况，为什么啊？

chybeyond · 发表于 2015-1-28 19:48

jerry0204103 发表于 2015-1-28 17:20
??? Attempt to execute SCRIPT STFT as a function.

Error in ==> STFTcs at 6

检查一下是不是有脚本文件命名为STFT，和函数STFT冲突

jerry0204103 · 发表于 2015-1-28 21:12

chybeyond 发表于 2015-1-28 19:48
检查一下是不是有脚本文件命名为STFT，和函数STFT冲突

function [Spec] = STFTwww(Sig,N,WinLen);

% Calculating the Short Time Fourier Transform
%
%          Sig    : the signal to be analyzed
%       N    : the number of frequency bins (default : length(Sig)).
%       WinLen    : the length of window used to locate the signal in time.
%
% June 1, 2003
%

if (nargin < 1),
error('At least one parameter required!');
end;

Sig = real(Sig);
SigLen = length(Sig);

if (nargin < 3),
WinLen = SigLen / 4;
end

if (nargin < 2),
N = SigLen;
end

if(N > 512),
   N = 512;
end

WinLen = ceil(WinLen / 2) * 2;
WinFun = exp(log(0.005) * linspace(-1,1,WinLen)'.^2 );
WinFun = WinFun / norm(WinFun);

Lh = (WinLen - 1)/2;
Spec = zeros(N,SigLen) ;

wait = waitbar(0,'Under calculation, please wait...');
for iLoop = 1:SigLen,

waitbar(iLoop/SigLen,wait);
tau = -min([round(N/2)-1,Lh,iLoop-1]):min([round(N/2)-1,Lh,SigLen-iLoop]);
Index = floor(rem(N + tau,N) + 1);
temp = floor(iLoop + tau);
temp1 = floor(Lh+1+tau);
Spec(Index,iLoop) = Sig(temp) .* conj(WinFun(temp1));

end;

close(wait);

Spec = fft(Spec);
Spec = abs(Spec);
f = [0:N/2-1 -N/2:-1]'/N;
t = 1: SigLen;

fmax = 0.5;
fmin = 0.0;

clf
%=====================================================%
% Plot the result                                  %
%=====================================================%
set(gcf,'Position',[20 100 800 500]);
set(gcf,'Color','w');

axes('Position',[0.07 0.37 0.4 0.55]);
mesh(t,f,Spec);
axis([min(t) max(t) fmin fmax]);
ylabel('Frequency');
title('Short Time Fourier Transform');

%=======================================%
% Plot the signal in the local zone    %
%=======================================%
axes('Position',[0.07 0.08 0.4 0.2]);
plot(t,real(Sig),'r');
xlabel('Time');
ylabel('Amp');
axis([1 SigLen min(real(Sig)) max(real(Sig))]);
title('Orignal signal');
grid

%===================================================%
% plot the FFT of the signal                      %
%===================================================%
sptr = fft(Sig);                                     %sptr is the FFT spectrum of the signal
sptrLen = length(sptr);
sptr = sptr(1:round(sptrLen/2)+1);
f = linspace(0,0.5,round(sptrLen/2)+1);
axes('Position',[0.55 0.37 0.4 0.55]);
plot(f,abs(sptr),'b');
axis([min(f) max(f) 0 max(abs(sptr))]);
title('Discrete Fourier Transform');
xlabel('Frequency');
ylabel('Amp');
grid

测试程序

SampFreq =30;
t=0:1/SampFreq:40;
sig = sin(12*pi*t);
sig(1:end/2) = sig(1:end/2) + sin(6*pi*t(1:end/2));
sig(end/2+1:end) = sig(end/2+1:end) + sin(18*pi*t(end/2+1:end));
STFTwww('sig',512,512);
这样改还有毛病，为啥？

chybeyond · 发表于 2015-1-28 21:36

jerry0204103 发表于 2015-1-28 21:12
function [Spec] = STFTwww(Sig,N,WinLen);

% Calculating the Short Time Fourier Transform

报什么错？STFTwww(sig',512,512);不STFTwww('sig',512,512);

jerry0204103 · 发表于 2015-1-28 21:46

chybeyond 发表于 2015-1-28 21:36
报什么错？STFTwww(sig',512,512);不STFTwww('sig',512,512);

这回对了，你给我的是正确的回答！十分感谢!

woshiqiao · 发表于 2015-5-25 21:45

账号		自动登录	找回密码
密码			我要加入

[小波] short time fourier transform

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

评分

测试程序

评分