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

short time fourier transform

function = 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():min();
    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 = '/N;
t = 1: SigLen;

fmax = 0.5;
fmin = 0.0;

clf
%=====================================================%
% Plot the result                                     %
%=====================================================%
set(gcf,'Position',);          
set(gcf,'Color','w');                                          

axes('Position',);             
mesh(t,f,Spec);
axis();
ylabel('Frequency');
title('Short Time Fourier Transform');

%=======================================%
% Plot the signal in the local zone   %
%=======================================%
axes('Position',);          
plot(t,real(Sig),'r');
xlabel('Time');
ylabel('Amp');
axis();
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',);
plot(f,abs(sptr),'b');
axis();
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 = 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():min();
    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 = '/N;
t = 1: SigLen;

fmax = 0.5;
fmin = 0.0;

clf
%=====================================================%
% Plot the result                                     %
%=====================================================%
set(gcf,'Position',);            
set(gcf,'Color','w');                                          

axes('Position',);               
mesh(t,f,Spec);
axis();
ylabel('Frequency');
title('Short Time Fourier Transform');

%=======================================%
% Plot the signal in the local zone   %
%=======================================%
axes('Position',);         
plot(t,real(Sig),'r');
xlabel('Time');
ylabel('Amp');
axis();
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',);
plot(f,abs(sptr),'b');
axis();
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 = 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

{:{39}:}{:{39}:}
页: [1]
查看完整版本: short time fourier transform