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

SlaveOpen

已于 2022-03-08 10:55:54 修改

阅读量4.5k

点赞数 8

文章标签： matlab 算法开发语言

于 2022-03-07 10:46:27 首次发布

本文链接：https://blog.csdn.net/SlaveOpen/article/details/123302158

版权

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

互相关算法原理
开发板实现互相关算法
DSP结果与MATLAB结果比较
LFM线性调频信号的互相关结果比较
数据及所有程序工程算法资料下载

互相关算法原理

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

互相关时延估计是一种最基本的时延估计算法，对接收到的两个声源信号进行互相关运算，互相关函数最大值对应两个声源信号到达麦克风的时间差τ。设x_o (n)和x_i (n)分别为FM参考传感器和麦克风i（i取1,2,3,4）接收到的信号，则他们之间的互相关为：
公式1互相关
对互相关公式1做FFT变换得：
公式2FFT变换
再对FFT变换公式2做IFFT变换得：
公式3IFFT变换
综上所述，得出互相关算法的流程图如下：
图2 互相关算法流程图
按照上述原理对音频信号做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

仿真结果如下：

matlab互相关算法求解输出结果

开发板实现互相关算法

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

由DFT和IDFT的公式可知：
公式4 DFT和公式5 IDFT
将公式4和公式5写成矩阵形式如下：
公式6 DFT矩阵形式和公式7 IDFT矩阵形式
公式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]);

FFT结果比较
IFFT结果比较
互相关结果比较
DSP互相关结果和MATLAB互相关结果作差比较

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)