数字锁相放大器（DLIA）基本原理与Matlab仿真-编程知识

数字锁相放大器（DLIA）基本原理与Matlab仿真

本文介绍数字锁相放大器（DLIA）基本原理与Matlab仿真。

1.基本原理

数字锁相放大器（DLIA）原理框图如下图。

其核心部分为FPGA/DSP内部的相关运算及由正弦参考序列，D/A转换器，低通滤波器构成的DDS，DDS相关内容可以参考我以前的博文，这里对信号相关运算的原理作简要介绍。

设输入信号为携带噪声的信号 $x(t)=A_{I}\cdot sin(2\pi ft+\theta )+n(t)$ ，输入信号频率为 $f$ ，采样频率为 $f_{s}$ ，令 $f_{s}=Nf(N\geq 3)$ ，采样间隔 $\tau =\frac{1}{Nf}$ ，对信号采样 $q$ 个周期，总采样点数为 $M=qN$ ， $n(t)$ 为噪声信号。参考信号分别为正弦序列和余弦序列（采样点数和输入信号一致），由于输入信号与参考信号之间有相关性，与噪声之间几乎无关，故在下面的相关运算过程中，噪声项可忽略。对输入信号进行采样可得离散信号序列：

$x(k)=A_{I}\cdot sin(2\pi fk\tau +\theta )=A_{I}\cdot sin(\frac{2k\pi}{N}+\theta) \, \, \, \, k=0,1,2... M-1$

参考信号可以由FPGA/DSP内部生成（采用表格法），这里分别生成正弦参考序列（ $r_{s}(k)$ ）和余弦参考序列（ $r_{c}(k)$ ）：

$r_{s}(k)=A_{R}\cdot sin(\tfrac{2k\pi}{N})\, \, \, \, k=0,1,2... M-1$

$r_{c}(k)=A_{R}\cdot cos(\tfrac{2k\pi}{N})\, \, \, \, k=0,1,2... M-1$

这里注意输入信号的频率和参考信号的频率要相同。

$x(k)$ 与 $r_{s}(k)$ 的互相关 $R_{xr,s}$ ， $x(k)$ 与 $r_{c}(k)$ 的互相关 $R_{xr,c}$ ，即同相输出和正交输出的互相关信号为：

$R_{xr,s}=\frac{1}{M-i}\sum_{k=0}^{M-i-1}x(k)\cdot r_{s}(k+i)=\frac{A_{I}A_{R}}{2}cos(2\pi fi+\theta)\, \, \, \, i=0,1,2...M-1$

$R_{xr,c}=\frac{1}{M-i}\sum_{k=0}^{M-i-1}x(k)\cdot r_{c}(k+i)=\frac{A_{I}A_{R}}{2}sin(2\pi fi+\theta)\, \, \, \, i=0,1,2...M-1$

其中， $R_{xr,s}$ ， $R_{xr,c}$ 采用的是无偏估计的相关运算。

从而可求得：

$A_{I}(i)=\frac{2{\sqrt{R_{xr,s}^{2}(i)+R_{xr,c}^{2}(i)}}}{A_{R}}\, \, \, \, i=0,1,2...M-1$

$\theta (0)=arctan(\frac{R_{xr,c}(0)}{R_{xr,s}(0)})$

2.Matlab仿真

这里，我们以输入信号fm=1000Hz，信噪比为-10dB（即噪声比信号强10倍），fs=100000Hz，为例，对上述原理进行仿真，Matlab代码如下：

fm=1000;
fs=100000;
N=fs/fm;
k=0:1:1000;
theta=pi/8;
x=sin(2*k*pi/N + theta);%原始信号
xn=awgn(x,-10,0);%对信号加噪
subplot(9,1,1);
plot(k,x);
title('x');
subplot(9,1,2);
plot(k,xn);
title('xn');
rs=sin(2*k*pi/N);%正弦信号
subplot(9,1,3);
plot(k,rs);
title('rs');
rc=cos(2*k*pi/N);%余弦信号
subplot(9,1,4);
plot(k,rc);
title('rc');rxs=xcorr(xn,rs,'unbiased',500);%与正弦信号作互相关
subplot(9,1,5);
plot(k,rxs);
title('rxs');
rxc=xcorr(xn,rc,'unbiased',500);%与余弦信号作互相关
subplot(9,1,6);
plot(k,rxc);
title('rxc');n = 10;
Wn = 0.3;
[b, a] = butter(n, Wn);rxs_filtered = filter(b, a, rxs);%对相关信号结果滤波
subplot(9,1,7);
plot(k,rxs_filtered);
title('rxs filtered');
rxc_filtered = filter(b, a, rxc);%对相关信号结果滤波
subplot(9,1,8);
plot(k,rxc_filtered);
title('rxc filtered');amp=2*sqrt(rxs_filtered.^2+rxc_filtered.^2);%计算幅值
subplot(9,1,9);
plot(k,amp);
title('amp');
p=atan(rxc_filtered(1)/rxs_filtered(1));%计算相位

其中，最核心的部分为相关运算（这里采用的是“unbiased”即无偏估计），相关运算计算可参考我以前的博文，这里不做过多叙述。仿真结果如下图。