Zoom FFT的matlab程序

songzy41 · 发表于 2006-8-31 18:29

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

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

x

ZoonFFT取自于Mathworks User’s Community中的程序，由Tom Irvine编制的。它的使用方法:在Matlab Command Window中输入
>>zoomfft
将显示出要求输入数据的方式（文件或已在Matlab中），若是文件要求输入文件名。如果以文件输入，文件中必须包含有时间和信号两部分：:
zoomFFT - version 1.1, April 20, 2005
  by Tom Irvine
  Email:  tomirvine@aol.com

This program calculates the non-destructive zoom
Fast Fourier transform of a time history.

The input file must be time(sec) and amplitude(units)
The format is free, but no header lines allowed.
Please enter the input filename.

Select file input method
1=external ASCII file
2=file preloaded into Matlab

Enter the input filename

再要求输入你希望细化的（中心）频率（即细化以该频率为中心）
Enter frequency (Hz) of interest

还要求输入细化时扩大的比例
Choose Zoom Factor
1 =  1:1
2 =  2:1
3 =  4:1
4 =  8:1
5 = 16:1

输入完成后便开始运行了，计算后给出以细化频率为中心的谱图。

ZoonFFT程序代码如下：

disp(' zoomFFT - version 1.1, April 20, 2005 ');
disp(' by Tom Irvine ');
disp(' Email: tomirvine@aol.com ');
%
disp(' ');
disp(' This program calculates the non-destructive zoom ');
disp(' Fast Fourier transform of a time history. ');
%
disp(' ');
disp(' The input file must be time(sec) and amplitude(units) ');
disp(' The format is free, but no header lines allowed.');
disp(' Please enter the input filename. ');
%
% Add phase compensation in future version
%
% Reference: Randall, Frequency Analysis, Bruel & Kjaer, page 170.
%
% 262144
% 131072
% 65536
% 32768
% 16384
%
clear THM;
clear amp;
clear ampf;
clear tim;
clear FF;
clear n;
clear N;
clear MM;
clear nz;
clear Y;
%
MAX = 262144*2;
MAXP1 = 131072*2;
%
inv=0;
mf=0.;
%
tp=2.*pi;
%
iflag=0;
%
disp(' ')
disp(' Select file input method ');
disp(' 1=external ASCII file ');
disp(' 2=file preloaded into Matlab ');
file_choice = input('');
%
if(file_choice==1)
disp(' Enter the input filename ');
filename = input(' ','s');
fid = fopen(filename,'r');
THM = fscanf(fid,'%g %g',[2 inf]);
THM=THM';
else
THM = input(' Enter the matrix name: ');
end
%
%N=4096;
%
amp=THM(:,2);
tim=THM(:,1);
n = length(amp);
%
for(i=1:18)
if( 2^i > n )
break;
end
N=2^i;
end
%
out4 = sprintf(' time history length = %d ',n);
disp(out4)
disp(' ');
out5 = sprintf(' samples used for FFT = %d',N);
disp(out5)
nsegment=N;
%
%disp(' mean values ')
%
mu=mean(amp);
sd=std(amp);
mx=max(amp);
mi=min(amp);
rms=sqrt(sd^2+mu^2);
tmx=max(tim);
tmi=min(tim);
dt=(tmx-tmi)/(n-1);
%
sr=1/dt;
df=1./(N*dt);
%
disp(' ')
out5 = sprintf(' dt=%12.4g sec sr=%12.4g Hz',dt,sr);
disp(out5)
disp(' ')
out5 = sprintf(' df=%12.4g Hz',df);
disp(out5)
%
%*** choose frequency
%
disp(' ');
disp(' Enter frequency (Hz) of interest ');
fin = input(' ');
%
octave = (2.^(1./3.));
flow = fin/octave;
fhigh = fin*octave;
%
%*** Choose Zoom Factor
%
disp(' ');
disp(' Choose Zoom Factor ');
disp(' 1 = 1:1 ');
disp(' 2 = 2:1 ');
disp(' 3 = 4:1 ');
disp(' 4 = 8:1 ');
disp(' 5 = 16:1 ');
%
iz=input(' ');
%
if(iz>5)
iz=5;
end
%
if(iz==1)
nz=1;
end
if(iz==2)
nz=2;
end
if(iz==3)
nz=4;
end
if(iz==4)
nz=8;
end
if(iz==5)
nz=16;
end
%
%***** Choose filter option (next version)
%
%disp(' ');
%disp(' Bandpass filter the data prior to zoom FFT (recommended) ?');
%disp(' 1=yes 2=no ');
%
%i_filter_choice=input(' ');
%
%if(i_filter_choice ==1)
%
% put in filter
%
%end
%
% FFT
%
disp(' ')
disp(' begin FFT ')
disp(' ')
%
clear Y;
clear complex_FFT;
Y=zeros(N,1);
complex_FFT=zeros(N,3);
%
ia=1;
MM=N;
N=fix(N/nz);
ja=1;
for( ikj = 1:nz)
nk = ikj;
jb=ja+MM-1;
ampf=zeros(1,MM);
%
ampf=amp(nk:nz:MM);
%
Y(ja:jb)=fft(ampf,MM);
ja=jb+1;
end
%
NT=length(Y);
FF=linspace(0,df*NT,NT);
complex_FFT = zeros(NT,3);
complex_FFT(:,1)=FF';
complex_FFT(:,2)=real(Y);
complex_FFT(:,3)=imag(Y);
%
nntt = fix(nsegment/nz);
%
% accumulator
%
disp(' accumulator ');
%
maxf = max(FF);
%
if( NT > MAX )
disp(' Error: too many data points. ');
disp(' Press any key to exit. ');
end
%
nmax=NT;
%
if( maxf < flow )
disp('\n Frequency error. Greater bandwidth needed. \n');
disp('\n Press any key to continue. \n');
input(' ');
end
%
dfz=df/nz;
nm= fix(NT/nz);
%
if(nz==1)
%
cr = complex_FFT(:,2);
ci = complex_FFT(:,3);
z=sqrt(cr.*cr+ci.*ci);
freq=FF;
%
else
%
for(i=1:nm)
%
cr=0.;
ci=0.;
for(j=0:(nz-1))
cr=cr + complex_FFT(i+j*nm,2) ;
ci=ci + complex_FFT(i+j*nm,3) ;
end
%
z(i)=sqrt(cr^2+ci^2);
freq(i)=i*dfz;
end
end
%
% Output
%
ny=fix(length(z)/2);
zoom=2.*z(1:ny);
freqz=freq(1:ny);
zoom(1)=zoom(1)/2.;
plot(freqz,zoom);
xlabel('Frequency (Hz)');
ylabel('Magnitude');
out5 = sprintf(' ZOOM FFT %d:1 ',nz);
title(out5);
grid;
fmin=fin/2.^(1./6.);
fmax=fin*2.^(1./6.);
ymin=0.;
ymax=max(zoom)*1.2;
axis([fmin,fmax,ymin,ymax]);

复制代码

[ 本帖最后由 suffer 于 2006-11-8 16:32 编辑 ]

mofei · 发表于 2006-9-5 20:09

我是新手，请教个问题：通过FFT可求出所分析信号的实部、虚部。也可以据此求出模。那么在频谱分析中用的是实部还是模。谢谢了！

shirley329 · 发表于 2006-11-7 11:42

取模
abs(shiftfft(fft(x))

realhappy · 发表于 2006-11-7 15:30

f=fl:df:fl+(N/2-1)*df;是否有错，应该是f=fl:df:fh+(N/2-1)*df;吧

songzy41 · 发表于 2006-11-7 20:16

奇怪，不是我上传的文件，而是
http://forum.vibunion.com/forum/ ... p;highlight=zoomfft
中的程序？？现重新把MATHWORKS的文件传上了。

suffer · 发表于 2006-11-8 16:33

原帖由 songzy41 于 2006-11-7 20:16 发表
奇怪，不是我上传的文件，而是
http://forum.vibunion.com/forum/ ... p;highlight=zoomfft
中的程序？？现重新把MATHWORKS的文件传上了。

呵呵，这就不清楚了！

realhappy · 发表于 2006-11-9 16:05

原帖由 realhappy 于 2006-11-7 15:30 发表
f=fl:df:fl+(N/2-1)*df;是否有错，应该是f=fl:df:fh+(N/2-1)*df;吧

怎么没有人回答这个问题呢，f是细化后的频率吗？

songzy41 · 发表于 2006-11-9 19:47

realhappy主任，你在4楼提到的问题，实际上不是我上传程序中的问题。我在5楼中指出，不知怎么搞的我上传的程序变成
http://forum.vibunion.com/forum/vi ... p;highlight=zoomfft
帖子上的程序。你提的问题是针对该帖子的。

jiangpq88 · 发表于 2006-11-29 14:10

好，很好

HolySaint · 发表于 2006-12-7 19:46

是ZoomFFT
文件名有个地方写错了

suffer · 发表于 2006-12-9 09:41

原帖由 HolySaint 于 2006-12-7 19:46 发表
是ZoomFFT
文件名有个地方写错了

这个不影响程序的运行

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

[共享资源] Zoom FFT的matlab程序

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

评分

本帖被以下淘专辑推荐: