gusucode.com > GPS相关检测源码程序 > GPS相关检测源码程序/fdcorr/fdcorr.m
function [fdout,freq_axis]= fdcorr(a,b,frange) %FDCORR Joint frequency and delay correlation % Performs an efficient circular cross correlation between % two vectors over all possible frequency offsets (Dopplers) % between -pi and pi where 2pi is the sampling rate. % %[FDOUT, FREQ_AXIS] = FDCORR(A,B) % % A: vector representing received signal sampled at FS = 2pi % B: vector representing correlation code sampled at FS = 2pi % % FDOUT: 2D Correlation surface vs delay and frequency offset % FREQ_AXIS: Frequency values aligned to frequency axis (-pi to pi) % %[FDOUT, FREQ_AXIS] = FDCORR(A,B,FRANGE) % % Limits frequency range as defined in FRANGE. Typically % frequency offsets are much smaller than the Nyquist bandwidth, so % using a reduced range for FRANGE from -pi to pi will be more % efficient. % % FRANGE: Two element vector from lowest to highest frequency to % search over within a frequency range from -pi to pi. % Example [-pi/10 pi/10] % %FDCORR(A,B,<FRANGE>) with no output arguments plots the surface % correlation magnitude vs frequency and delay offset in the % current figure window. FRANGE optional. % % % Dan Boschen 11/22/2009 %========================================================== % VERSION HISTORY % % 11/22/2009 Inital Release % % %========================================================== %========================================================== % APPLICATION NOTES: % % Vector "A" represents a received signal that has been encoded with % a correlation code (such as GPS and CDMA applications) or a training % sequence. It has an unknown frequency offset due to Doppler offsets % and frequency offsets in the receiver, and an unkown delay shift % relative to the start of the code sequence represented by vector "B". % FDCORR will efficiently find the correlation between the received % sequence and reference code over frequency and delay offsets. % Vectors A and B are interchangable. % % To detect the maximum correlation value and delay and freq location use: % [temp1, temp2]=max(fdout); % [maxcorr,freq_index]=max(temp1); % delay =temp2(freq_index)-1; % freq= freq_axis(freq_index); % % The basic operation uses the Cross-Correlation theorem: % xcorr= ifft(fft(a).*ifft(b)) % % Note the detected frequency offset is course and is the % frequency of the closest Doppler bin where each Doppler bin % is seperated by 2pi/N % %========================================================== % %error handling if nargin==3 if length(frange)~= 2 error('FRANGE must be a two element vector') elseif abs(frange)> pi error('FRANGE values must be between -pi and pi') end end %force column vectors: a=a(:); b=b(:); %============================================================ %create all doppler versions of a by shifting the fft: %============================================================ N=max(length(a),length(b)); fft_a= fft(a,N); %frequency axis for all doppler bins -pi to pi %index is the "fft shifted" freq axis. Method below faster than %using fftshift(0:N-1), but same result. freq_index=[-floor((N-1)/2):floor(N/2)]; freq_axis= freq_index*2*pi/N; %convert to -pi to pi range if nargin==3 %reduced doppler range freq_include=find(freq_axis>=frange(1) & freq_axis<=frange(2)); freq_axis=freq_axis(freq_include); freq_index=freq_index(freq_include); num_freqs=length(freq_index); else num_freqs=N; end %create Hankel matrix for all Dopplers (each column is the fft %of vector a at a different Doppler offset) fft_idx=(1:N)'; fft_dopplers = fft_idx(:,ones(num_freqs,1))+freq_index(ones(N,1),:); fft_dopplers= mod(fft_dopplers-1,N)+1; %Hankel subscripts fft_dopplers(:)= fft_a(fft_dopplers); %fill data %alternate method if solving for all Doppler bins (processing intensive): %fft_all_dopplers= hankel(fft_a,[fft_a(end); fft_a(2:end)]); %================================================================== %Circular Correlation using circ_corr= ifft(fft(a).*conj(fft(b)) %================================================================== fft_code_mtx= fft(b,N)*ones(1,num_freqs); fdout=(ifft(fft_code_mtx.*conj(fft_dopplers))); %=============================================================== %plot result %=============================================================== if nargout ==0 surf(freq_axis,1:N,(abs(fdout))); shading interp; xlabel('Normalized Frequency (xpi rad/sample)'); ylabel('Delay Index'); end