互相关算法原理

为了实现声源定位，需要求出阵列之间的时延信息，常用的方式就是互相关算法。因为定位系统需要机动灵活的特性的特性，互相关算法的实现不能完全依赖于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:4096fprintf(fid,'%f\r\n',y(1,jj));
end
fclose(fid)

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

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

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

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

DSP库互相关算法实现与MATLAB互相关算法比较相关推荐

DFP算法求极值点matlab,DFP算法及Matlab程序
作业二用DFP 算法求解1212221422)(m in x x x x x x f --+=,取()T x 110=,??? ? ??=10010H . 一.求解: T T T g H p g x ...
MATLAB实战系列（十一）-多种群遗传算法的函数优化算法（附MATLAB代码）
前言: 本篇博文参考,智能优化算法书籍<MATLAB智能算法30个案例分析(第2版)>,今天要与大家分享的智能算法是多种群遗传算法. 本地MATLAB环境部署因为后面要介绍的多种群遗传算 ...
粒子群算法(PSO)以及Matlab实现
粒子群算法(PSO)以及Matlab实现算法背景粒子群优化算法(PSO)是一种进化计算技术(evolutionary computation),1995 年由Eberhart 博士和kennedy ...
stm32使用dsp库，结合Matlab进行FIR滤波器设计
首先我们打开Matlab.在命令串口输入fdatool,按回车. Response Type :这里可以设置滤波器类型 lowpass(低通),highpass(高通),banpass(带通),ban ...
SVPWM仿真和基于DSP28335的PIL(处理器在环) 仿真模型（将matlab仿真算法生成代码在DSP中在线运行返回数据给Matlab）验证算法可行性和实时性
SVPWM仿真和基于DSP28335的PIL(处理器在环) 仿真模型(将matlab仿真算法生成代码在DSP中在线运行返回数据给Matlab)验证算法可行性和实时性. 对于数字信号处理很有用. ID: ...
matlab 图像拼接算法,MATLAB图像拼接算法及实现
MATLAB图像拼接算法及实现图像拼接算法及实现(一)论文关键词:图像拼接图像配准图像融合全景图论文摘要:图像拼接(image mosaic)技术是将一组相互间重叠部分的图像序列进行空间匹配对 ...
蚁群算法优化神经网络matlab源程序,粒子群优化神经网络的程序大集合
粒子群程序集合 866867259psobp psobp.m pso(粒子群算法)优化神经网络粒子群算法(PSO)应用于神经网络优化[matlab] PSOt A Particle Swarm Op ...
使用FDDB人脸样本检测库，测试自己的人脸检测算法性能并生成ROC曲线。
一,说明及环境网上有关FDDB人脸检测库的使用以及ROC文件生成的文章太少,并且都无法检测opencv中自带的人脸检测算法.最近工作的原因,需要用到FDDB库检测我们自己的人脸检测算法性能.所以认 ...
【故障检测问题】基于matlab免疫算法求解故障检测问题【含Matlab源码 196期】
一.获取代码方式获取代码方式1: 完整代码已上传我的资源:[故障检测问题]基于matlab免疫算法求解故障检测问题[含Matlab源码 196期] 获取代码方式2: 通过订阅紫极神光博客付费专栏,凭 ...

DSP库互相关算法实现与MATLAB互相关算法比较