zoomfft.m

matlab中的一个zoomfft例程

function Sigzoomffted = zoomfft(fid,fin,iz)

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;
    THM = fid;
    

%
amp=THM(:,2);   %取第二列作为信号的幅值
tim=THM(:,1);   %取第一列作为信号的时间
n = length(amp);
%取小于且与n最接近的2的整数次幂
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);    %标准偏差S
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         
%        

%确定频率范围
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 ');
%

%
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Y0 = Y(1:MM)+Y(MM+1:2*MM)+Y(2*MM+1:3*MM)+Y(3*MM+1:4*MM);
%Y0 = Y0./4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NT=length(Y);

out9 = sprintf(' tht length of NT is NT = %d ',NT);
disp(out9);

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)=cr+i*ci;%
        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.;
%Sigzoomffted = zoom;
Sigzoomffted = Y0;




%plot(freqz,zoom);
%xlabel('Frequency (Hz)');
%ylabel('Magnitude');
%out5 = sprintf(' ZOOM FFT  %d:1 ',nz);
%title(out5);
%grid;
%fmin=fin/30.^(1./6.);
%fmax=fin*30.^(1./6.);
%ymin=0.;
%ymax=max(zoom)*1.2;
%axis([fmin,fmax,ymin,ymax]);