short time fourier transform

pengzk

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

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

多谢共享自己的成果

songzy41

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

ghb063

非常感谢

jerry0204103

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

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

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

chybeyond

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

Error in ==> STFTcs at 6

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

jerry0204103

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

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

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

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

woshiqiao