互相关算法原理

为了实现声源定位，需要求出阵列之间的时延信息，常用的方式就是互相关算法。因为定位系统需要机动灵活的特性的特性，互相关算法的实现不能完全依赖于PC端软件，所以本文尝试了基于开发板DSP库的互相关算法实现。

互相关时延估计是一种最基本的时延估计算法，对接收到的两个声源信号进行互相关运算，互相关函数最大值对应两个声源信号到达麦克风的时间差τ。设x_o (n)和x_i (n)分别为FM参考传感器和麦克风i（i取1,2,3,4）接收到的信号，则他们之间的互相关为：

对互相关公式1做FFT变换得：

再对FFT变换公式2做IFFT变换得：

综上所述，得出互相关算法的流程图如下：

按照上述原理对音频信号做matlab的仿真，仿真条件如下：
假设已知信号为起始频率为250Hz，终止频率为2000Hz的Chirp信号，接收信号为麦克风利用ADC采集到的在幅度上有一定衰减的Chirp信号。

%%% 互相关法求时延 %%%
N=2048;             % 长度
Fs=10000;           % 采样频率
n=0:N-1;
t=n/Fs;             % 时间序列
%% Chirp信号
x = chirp(t,250,0.2047,2000);
x = [x,zeros(1,2*N)];
                    % 要保证时延在两倍信号长度之内，不然就调整2*N的2
%% 接收到的Chirp信号
ADC1 = 0.78 * chirp(t,250,0.2047,2000);
                    % 接收后的信息假设存在幅度衰减
tao = 0.15;         % 延时
Ndelay = fix(tao*Fs);
                    % 换算成延时点数
ADC1 = [zeros(1,Ndelay),ADC1,zeros(1,length(x)-Ndelay-length(ADC1))];%,zeros(1,2048)];
%% 互相关函数
[c,lags]=xcorr(x,ADC1);
                   % 互相关
subplot(211);
xaxis = 1:1:length(x);
plot(xaxis,x,'r');
hold on;
plot(xaxis,ADC1,':');
legend('x信号', 'ADC1信号');
xlabel('时间/s');ylabel('x(t) ADC1(t)');
title('原始信号');grid on;
hold off
subplot(212);
plot(lags/Fs,c,'r');
title('相关函数');
xlabel('时间(s)');ylabel('Rxadc1');
grid on
[Am,Lm]=max(c);
d = Lm - (length(c)+1)/2;
phy=(2*10*d*180/Fs);
phy = rem(phy,360)% 取余360
phy =
  -180
Delay=d/Fs
Delay =
   -0.1500

仿真结果如下：

开发板实现互相关算法

离散傅里叶变换（DFT）与离散傅里叶反变换（IDFT）

由DFT和IDFT的公式可知：

将公式4和公式5写成矩阵形式如下：

公式4和公式5中的W为旋转因子，只是DFT和IDFT的旋转因子W是共轭的，所以对于DFT和IDFT，都相当于进行了旋转因子不同的DFT运算，由欧拉公式可知：

DFT旋转因子的计算

void tw_gen(float *w, int n)
{
    int i, j, k;
    const double PI = 3.141592654;
    for (j = 1, k = 0; j <= n >> 2; j = j << 2) {
        for (i = 0; i < n >> 2; i += j) {
            w[k]     = (float) sinsp (2 * PI * i / n);
            w[k + 1] = (float) cossp (2 * PI * i / n);
            w[k + 2] = (float) sinsp (4 * PI * i / n);
            w[k + 3] = (float) cossp (4 * PI * i / n);
            w[k + 4] = (float) sinsp (6 * PI * i / n);
            w[k + 5] = (float) cossp (6 * PI * i / n);
            k += 6;
        }
    }
}

IDFT旋转因子的计算

根据欧拉公式，只需要对DFT的旋转因子做共轭变换即可。

void tw_geni(float *w, int n)
{
    int i, j, k;
    const double PI = 3.141592654;
    for (j = 1, k = 0; j <= n >> 2; j = j << 2) {
        for (i = 0; i < n >> 2; i += j) {
            w[k]     = (float) cossp (2 * PI * i / n);
            w[k + 1] = (float) -sinsp (2 * PI * i / n);
            w[k + 2] = (float) cossp (4 * PI * i / n);
            w[k + 3] = (float) -sinsp (4 * PI * i / n);
            w[k + 4] = (float) cossp (6 * PI * i / n);
            w[k + 5] = (float) -sinsp (6 * PI * i / n);
            k += 6;
        }
    }
}

FFT变换函数和IFFT变换函数

void FFT(float *Input, float *Cmo, int Tn)
{
    int i;
    int rad;
    if(Tn==16 || Tn==64 || Tn==256 || Tn==1024 || Tn==4096 || Tn==16384 || Tn==65536) {
        rad=4;
    } else if(Tn==8 || Tn==32 || Tn==128 || Tn==512 || Tn==2048 || Tn==8192 || Tn==32768 || Tn==131072) {
        rad=2;
    } else {
        return;
    }

    for(i=0;i<Tn;i++) {
        CFFT_In[2*i]=Input[i];
        CFFT_In[2*i+1]=0;
    }

    tw_gen(Cw,Tn);

    DSPF_sp_fftSPxSP(Tn,CFFT_In,Cw,CFFT_Out,brev,rad,0,Tn);

    for(i=0;i<Tn;i++) {
        if(i==0) {
            Cmo[i]=sqrtsp(CFFT_Out[2*i]*CFFT_Out[2*i]+CFFT_Out[2*i+1]*CFFT_Out[2*i+1])/Tn;
        } else {
            Cmo[i]=sqrtsp(CFFT_Out[2*i]*CFFT_Out[2*i]+CFFT_Out[2*i+1]*CFFT_Out[2*i+1])*2/Tn;
        }
    }
}

void IFFT(float *Input, float *Cmo, int Tn)
{
    int i;
    int rad;
    if(Tn==16 || Tn==64 || Tn==256 || Tn==1024 || Tn==4096 || Tn==16384 || Tn==65536) {
        rad=4;
    } else if(Tn==8 || Tn==32 || Tn==128 || Tn==512 || Tn==2048 || Tn==8192 || Tn==32768 || Tn==131072) {
        rad=2;
    } else {
        return;
    }

    for(i=0;i<Tn;i++) {
        CIFFT_In[2*i]=Input[i];//瀹為儴
        CIFFT_In[2*i+1]=0;//铏氶儴
    }

    tw_geni(CwI,Tn);

    DSPF_sp_fftSPxSP(Tn,CIFFT_In,CwI,CIFFT_Out,brev,rad,0,Tn);

    for(i=0;i<Tn*2;i++)
    {
        Cmo_IFFTOut[i]=CIFFT_Out[i]/Tn;
    }

 /*   for(i=0;i<Tn;i++)
    {
        Cmo_IFFTOut[2*i]= Cmo_IFFTOut[2*i];//瀹為儴
        Cmo_IFFTOut[2*i+1]=-Cmo_IFFTOut[2*i+1];//铏氶儴
    }*/
    for(i=0;i<Tn;i++)
    {
        if(i==0)
        {
            Cmo[i]=sqrtsp(Cmo_IFFTOut[2*i]*Cmo_IFFTOut[2*i]+Cmo_IFFTOut[2*i+1]*Cmo_IFFTOut[2*i+1]);
        }
        else
        {
            Cmo[i]=sqrtsp(Cmo_IFFTOut[2*i]*Cmo_IFFTOut[2*i]+Cmo_IFFTOut[2*i+1]*Cmo_IFFTOut[2*i+1]);
        }
    }
}

DSP结果与MATLAB结果比较

本文利用host端产生两个频率叠加正弦信号（4096点）：

  for(j = 0; j < PAYLOADSIZE; j++) {
        dataIn[j] = sin (2 * 3.1415 * 50 * j / (double) PAYLOADSIZE);
    }

DSP处理完成后发给host端，host端将结果存入文件中以进行后续处理。

close all;
clear;
clc;
time_txt = textread('C:\Users\lenovo\Desktop\TimeData');
fft_txt = textread('C:\Users\lenovo\Desktop\FFTData');
ifft_txt = textread('C:\Users\lenovo\Desktop\IFFTData');
rfft_txt = textread('C:\Users\lenovo\Desktop\RFFTData');

fft_local=2*(abs((fft(time_txt))))/4096;
ifft_local = abs((ifft(time_txt)));
LEN = 4096;
figure;
subplot(211);
plot(0:1:LEN-1,fft_txt(1:1:LEN));
title('DSP FFT')
subplot(212);
plot(0:1:LEN-1,fft_local(1:1:LEN));
title('MATLAB FFT')
figure;
subplot(211);
plot(0:1:LEN-1,ifft_txt(1:1:LEN));
title('DSP IFFT')
subplot(212);
plot(0:1:LEN-1,ifft_local(1:1:LEN));
title('MATLAB IFFT')


X1=(fft(time_txt));
X2=(fft(time_txt));
Sxy=X1.*conj(X2);
Bxy=fftshift(abs(ifft(Sxy)));
figure;
subplot(211);
plot((rfft_txt));
title('DSP hxg')
subplot(212);
plot((Bxy));
title('MATLAB hxg')
[max,location]=max(Bxy);
%d=location-N/2-1        
% d=location-4096
% Delay=d/Fs              %求得时间延迟
% Xxy=xcorr(time_txt,time_txt);
% plot(Xxy);
figure;
DSS=rfft_txt-Bxy;
plot(DSS);
axis([0 4096-1 -1 1]);

LFM线性调频信号的互相关结果比较

利用matlab生成线性调频信号，并将信号保存到文件中。

fs = 1000;
T = 5;
B = 10;
k = B / T;
n = round(T * fs);
t = linspace(0,T,n);
y = exp(1j*pi*k*t.^2);
figure;
plot(t,y);
xlabel('t/s');
ylabel('幅度');
X1=(fft(y));
X2=(fft(y));
Sxy=X1.*conj(X2);
Bxy=fftshift(abs(ifft(Sxy)));
figure;
plot((Bxy));
fid = fopen('C:\Users\lenovo\Desktop\lfmDATA','a');
for jj = 1:1:4096
  fprintf(fid,'%f\r\n',y(1,jj));  
end
fclose(fid)

生成后的线性调频信号时域波形如下：

通过上述算法对线性调频信号进行互相关运算得出结果如下：

数据及所有程序工程算法资料下载

点击数据及所有程序工程算法资料下载

function C=corr(image,template)
    [mi,ni]=size(image);
    [mt,nt]=size(template);
    C=zeros(mi-mt+1,ni-nt+1);
    for u=1:mi-mt+1 
        for v=1:ni-nt+1
            a=0;
            b=0;
            c=0;
            for i=1:mt
                for j=1:nt
                    a=a+image(i+u-1,j+v-1)*template(i,j);
                    b=b+image(i+u-1,j+v-1)*image(i+u-1,j+v-1);
                    c=c+template(i,j)*template(i,j);
                end
            end
            b=sqrt(b);
            c=sqrt(c);
            C(u,v)=a/(b*c);
        end
    end
end

S为实际图像，T为参考图像。

算法公式为：

1.归一化互相关匹配

https://blog.csdn.net/u013049912/article/details/85984238；

2.图像匹配—NCC算法，即归一化互相关匹配

https://blog.csdn.net/qq_32261049/article/details/78666654

3.图像处理之积分图应用三（基于NCC快速相似度匹配算法）

https://blog.csdn.net/fangyan90617/article/details/101269473

4.图像处理之基于NCC模板匹配识别

https://blog.csdn.net/weixin_34208185/article/details/89936381
5.图像匹配之归一化积相关灰度匹配

https://blog.csdn.net/cui134/article/details/23924303

6.双目立体匹配流程——归一化互相关（NCC）

https://blog.csdn.net/aaronmorgan/article/details/79121434