TDOA定位中的广义相关（Generalized Cross Correlation）算法仿真

最近在优化频谱监测系统中的TDOA定位精度。TDOA中最重要的问题无外乎是信号时差的提取。看了国内外的文献，一致首推“广义相关（Generalized Cross Correlation--GCC）算法”。

查文献，在matlab中使用模拟正弦波仿真测试了一下广义相关算法，没有达到预期效果，百思不得其解。浏览了一下网络上关于GCC的评论，国内外不少人都抱怨无法得到预期仿真结果。最后，想到广义相关算法公式中，影响计算结果的主要参数是相位信息，应该需要使用仿真调制信号来测试。使用Matlab仿真FM/FSK/PSK之后，效果果然不错。在C++程序中实现这个算法只是工作量的问题了。matlab仿真源代码展示在这里，供大家参考。注意，不同版本Matlab需要使用不同的调制函数。运行结果下次再贴上。

%%%%%%广义互相关算法仿真   Generalized Cross Correlation %%%%%%%%%%%

%%%%%% 【信号仿真】  %%%%%%%%%%%

        clear;

        ModType = 'fm';  %%信号类型可选‘fm’和‘fsk’和‘psk’三种类型

        deltaT = -100;   %信号之间的数据点延迟--可正负

        fftpoints=1024; %number of fft points

        M = 2; freqsep = 4096; SamplesPerSymbol = 8; Fs = freqsep * M; %最恰当采样率应为M和频率步长乘积  

        InterpFactor = 4; SymbolRate = Fs/SamplesPerSymbol             %调制信号的插值因子

        if strcmp(ModType, 'fm') == 1

        Xsig = Func_SigGen_AmSsbFm('fm', 2, SamplesPerSymbol, SymbolRate, 100, 1, 10)    

        end

        if strcmp(ModType, 'fsk') == 1

            sig = randint(4 * fftpoints/SamplesPerSymbol,1,M);             % Random signal               

            XsigSrc = fskmod(sig,M,freqsep,SamplesPerSymbol,Fs);            % Modulate：数据点数量为4倍fftpoints

            Xsig = interp(XsigSrc,InterpFactor);                           %信号插值   

        end

         if strcmp(ModType, 'psk') == 1

            sig = randint(4 * fftpoints/SamplesPerSymbol,1,M);             % Random signal               

            XsigSrc = pskmod(sig,M);            % Modulate：数据点数量为4倍fftpoints

            Xsig = interp(XsigSrc,InterpFactor);                           %信号插值   

        end

       %ly = length(Xsig);

        % Create an FFT plot.

%         freq = [-Fs/2 : Fs/ly : Fs/2 - Fs/ly];

%         XsigFft = 10*log10(fftshift(abs(fft(Xsig))));

%         plot(freq,XsigFft)% 

        cnt = 1;

        for k= fftpoints/2 + deltaT: fftpoints + fftpoints/2 + deltaT

            BsigIn(cnt) = Xsig(k);                          %转储B路信号

            cnt = cnt + 1;

        end

        

        cnt = 1;

        for k= fftpoints/2: fftpoints + fftpoints/2

            AsigIn(cnt) = Xsig(k);                          %转储A路信号

            cnt = cnt + 1;

        end

        

        

        

        t = 1:1:fftpoints;       

        %x = sin(2*pi*0.005*t)+sin(2*pi*0.012*t); %第一个信号

        AsigIn0 = awgn(AsigIn,0);    %给信道 增加 噪声

        %t = deltaT:1:fftpoints+deltaT-1;

        %y = sin(2*pi*0.005*t)+sin(2*pi*0.012*t); %第二个信号

        BsigIn0 = awgn(BsigIn,0);    %给信道 增加 噪声

        

        

        figure(1);    clf;    hold on;

        plot(t(1:fftpoints),real(AsigIn0(1:fftpoints)), '--ro','LineWidth',1,...

                'MarkerEdgeColor','r',...

                'MarkerFaceColor','r',...

                'MarkerSize',2)

        plot(t(1:fftpoints),real(BsigIn0(1:fftpoints)), '--go','LineWidth',1,...

                'MarkerEdgeColor','g',...

                'MarkerFaceColor','g',...

                'MarkerSize',2)

        title('Simulation Signa');    xlabel('Frequency (Hz)');    ylabel('|Y(f)|');

        hold off;

%%%%%% 【算法01】  %%%%%%%%%%%----------GCC的实现--基于 互相关功率谱的算法--可以运行，对采样率较敏感，但效果并不好

%         fs= InterpFactor * Fs;  %%采样率 过高，将导致GCC的计算错误！！！！！！！！0.2* 48000是合适的-->结果和采样率相关

%         L=128; %window length

%         %Let's evaluate the CPSD by means of the Welch's Method

%         Pxy=cpsd(AsigIn0,BsigIn0,hamming(L),L/2,fftpoints,fs,'twosided');

%         %Take care of fftpoints value considering the frame length!!!

%         GCC = ifftshift(real(ifft(Pxy./abs(Pxy)))); %GCC-PHAT evaluation

% 

%         figure(2);    clf;    hold on;

%                 plot(t(1:fftpoints) -fftpoints/2, GCC, '--bo','LineWidth',1,...

%                         'MarkerEdgeColor','r',...

%                         'MarkerFaceColor','r',...

%                         'MarkerSize',2)

%                 title('Simulation Signa');    xlabel('Frequency (Hz)');    ylabel('|Y(f)|');

%                 hold off;

%%%%%% 【算法02】%%%%%% %%%%%%%%%%%----------GCC的实现--基于FFTpinpu的算法--这个算法的性能较好！！！

%%%%%%这个算法需要使用信号来验证，随机信号不能用于验证！

        b1f=fft(AsigIn0);

        b2f=fft(BsigIn0);

        b2fc=conj(b2f);

        neuma=(b1f).*(b2fc);

        deno=abs(neuma)    % 

        %deno=abs((b1f).*(b2fc));    % 等效于abs(neuma)

        GPHAT=neuma./deno;

        GPHATi=ifft(GPHAT);

        [maxval ind]= max(abs(GPHATi));

        samp=ind;

        

        figure(3);    clf;    hold on;

        plot(abs(GPHATi), '--bo','LineWidth',1,...

                'MarkerEdgeColor','r',...

                'MarkerFaceColor','r',...

                'MarkerSize',2)

        title('Simulation Signa');    xlabel('Frequency (Hz)');    ylabel('|Y(f)|');

        hold off;