www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/zoomfft.m
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]);