1200字范文 > DSP库互相关算法实现与MATLAB互相关算法比较

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

时间：2023-04-02 22:40:55

互相关算法原理开发板实现互相关算法离散傅里叶变换（DFT）与离散傅里叶反变换（IDFT）DFT旋转因子的计算IDFT旋转因子的计算FFT变换函数和IFFT变换函数DSP结果与MATLAB结果比较LFM线性调频信号的互相关结果比较数据及所有程序工程算法资料下载

互相关算法原理

为了实现声源定位，需要求出阵列之间的时延信息，常用的方式就是互相关算法。因为定位系统需要机动灵活的特性的特性，互相关算法的实现不能完全依赖于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 offsubplot(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)% 取余360phy =-180Delay=d/FsDelay =-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)); endfclose(fid)

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

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