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
测试程序
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); 多谢共享自己的成果 感谢共享,不过很遗憾的是STFT只能处理实数信号,因为在STFT函数中有
Sig = real(Sig);
语句。如果能处理复数信号则更完善了。 非常感谢 ??? Attempt to execute SCRIPT STFT as a function.
Error in ==> STFTcs at 6
STFT('sig',512,512);
运行出现上边的情况,为什么啊? jerry0204103 发表于 2015-1-28 17:20
??? Attempt to execute SCRIPT STFT as a function.
Error in ==> STFTcs at 6
检查一下是不是有脚本文件命名为STFT,和函数STFT冲突 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);
这样改还有毛病,为啥? 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); chybeyond 发表于 2015-1-28 21:36
报什么错?STFTwww(sig',512,512);不STFTwww('sig',512,512);
这回对了,你给我的是正确的回答!十分感谢! {:{39}:}{:{39}:}
页:
[1]