Python模拟锁相放大器的工作过程

准备

  本篇笔记将使用python演示锁相放大器的基本原理。
  首先引入需要用到的package,使用%matplotlib widget可以产生交互式的图片。

1
2
3
4
5
6
import numpy as np
from scipy.fftpack import fft
from scipy import signal as sg
from matplotlib import pyplot as plt

%matplotlib widget

绘图

  在本文中将会产生许多的图片,为了方便绘图编写一个展示信号的绘图函数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
def draw(x, y, title='', xlable='time(s)', ylable='amplitude', L=-3.8, H=3.8):
'''
draw

Parameters
----------
x : numpy array
abscissa
y : numpy array
ordinate

Returns
-------
'''
# canvas settings
plt.close()
if L!=H:
plt.ylim(L, H)
plt.title(title)
plt.xlabel(xlable)
plt.ylabel(ylable)

plt.plot(x, y)
plt.show()

正弦波

  在本实验中将采用正弦波作为输入信号。另外,参考信号、噪声信号等都需要用到正弦信号,所以第一步将编写一个产生给定幅度、频率和相位的正弦信号的函数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
def g_sin(t=1, amp=1, f0=10, fs=500, phi=0):
'''
generate sinusoid

Parameters
----------
t : float
time length of the generated sequence
amp : float
amplitude
f0 : float
frequency of sinusoid in Hz
fs : float
sampling rate per second
phi : float
initial phase in deg

Returns
-------
anonymous : list[numpy array, numpy array]
[abscissa, a sinusoid signal]
'''
T = 1/fs
N = t/T
x = 2*np.pi*np.arange(N)*T
return [np.arange(0, t, T), amp*np.sin(f0*x+phi*np.pi/180)]


x, y = g_sin()
draw(x, y)

方波

  参考信号可以使用正弦波、方波等。在本实验中将使用方波,而方波的傅里叶级数为:

f(x)=4πn=0sin[(2n+1)ωrt]2n+1f(x)=\frac{4}{\pi}\sum\limits_{n=0}^{\infty}\frac{sin[(2n+1)\omega_rt]}{2n+1}

  所以我们可以利用上面的正弦波函数,产生一个近似的方波,方波的阶数(即K)越大近似效果越好,K=50时就有很好的效果。 (PS:此处方波均指占空比为50%,正负对称的方波)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
def g_square_wave(t=1, amp=1, f0=10, fs=500, K=50):
'''
generate square wave

Parameters
----------
t : float
time length of the generated sequence
amp : float
amplitude
f0 : float
frequency of sinusoid in Hz
fs : float
sampling rate per second
K : int
order of fourier series

Returns
-------
anonymous : list[numpy array, numpy array]
[abscissa, a square wave]
'''
T = 1/fs
N = t/T
x = 2*np.pi*np.arange(N)*T

y = np.zeros(len(x))
k = 2*np.arange(1, K)-1
for i in range(len(x)):
y[i] = amp*np.sum(np.sin(k*f0*x[i])/k)
return [np.arange(0, t, T), 4*y/np.pi]


x, y = g_square_wave()
draw(x, y)

  使用上述两个函数,可以叠加产生带有噪声的信号。

1
2
3
4
x, y1 = g_sin(amp=0.2, f0=30)
x, y2 = g_sin(amp=0.35, f0=3, phi=45)
x, y3 = g_square_wave(amp=0.1, f0=10)
draw(x, y1+y2+y3)

频域分析

  这些噪声在时域中难以分辨。不过信号与噪声不同,往往具有明显的频率特征,而噪声一般是与频率无关的,或者是在特定频率范围内的,所以在频域内可以很好的分辨出输入信号的各种成分。
  使用快速傅里叶变换(FFT)从时域转化到频域,并且对振幅谱进行取半归一化。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def FFT(s, fs):
'''
fast fourier transform

Parameters
----------
s : numpy array
signal
fs : float
sampling rate per second

Returns
-------
'''
X = fft(s)
mX = np.abs(X) # magnitude
pX = np.angle(X) # phase

draw(np.arange(int(fs/2)), mX[range(int(fs/2))]/fs,
'Input signal in frequency domain', 'frequency(Hz)', 'amplitude', 0, 0)


FFT(y1+y2+y3, 500)

高斯白噪声

  上面的信号添加的噪声还是正弦信号,而现实中比较常见的噪声是白噪声。下面的函数可以根据输入信号产生一个特定信噪比的白噪声。
  信噪比的定义是:

SNR(dB)=10lgPsignalPnoise=20lgAsignalAnoiseSNR(dB)=10lg\frac{P_{signal}}{P_{noise}}=20lg\frac{A_{signal}}{A_{noise}}

  在本实验中,功率PP可以用功W=Pdt=CA2dtW=\int Pdt=C\int A^2dt代替,因为信号和噪声的积分时间是相同的。而信号的波形是采样得到的,所以WW可以用xi2\sum x_i^2代替,于是有:

SNR=10lgsi2ni2ni2=si210SNR/10SNR=10lg\frac{\sum s_i^2}{\sum n_i^2} \Longrightarrow \sum n_i^2= \frac{\sum s_i^2}{10^{SNR/10}}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
def wgn(x, snr):
'''
generate Gauss white noise

Parameters
----------
x : numpy array
signal
snr : float
signal-to-noise ratio

Returns
-------
anonymous : numpy array
Gauss white noise
'''
snr = 10**(snr/10.0)
xpower = np.sum(x**2)/len(x)
npower = xpower / snr
return np.random.randn(len(x)) * np.sqrt(npower)


noise = wgn(np.sin(np.arange(0, 1000000)*0.1), 6)
plt.close()
plt.subplot(211)
plt.hist(noise, bins=100)
plt.subplot(212)
plt.psd(noise)
plt.show()

模拟

  下面开始正式的模拟过程:

  1. 获得输入信号,输入信号由10Hz的目标信号和其他频率的正弦噪声以及白噪声组成

  A. 增加正弦部分

1
2
3
4
5
x, y1  = g_sin(amp=0.1)
x, y2 = g_sin(amp=0.5, f0=50)
x, y3 = g_sin(amp=0.8, f0=3)
x, y4 = g_sin(amp=0.3, f0=200)
signal = y1+y2+y3+y4

  B. 增加白噪声部分

1
2
3
noise = wgn(signal, 10)
sn = signal+noise
draw(x, sn)


  C. 查看输入信号的频域信息

1
FFT(sn, 500)

  1. 获得参考信号,方波频率与目标信号相同,为10Hz
1
2
3
x, y   = g_square_wave(K=100)
square = y
draw(x, square)

  1. 前置放大×2倍
1
2
sn = sn*2
draw(x, sn)

  1. 输入信号进行带通滤波,除去/抑制部分噪声,增强锁相放大器的性能(动态储备和输出动态范围)
1
2
3
b, a = sg.butter(1, [0.02, 0.12], 'bandpass') # 5Hz - 30Hz
snq = sg.filtfilt(b, a, sn)
draw(x, snq)

  1. 调谐放大输入信号×5,增大信噪比
1
2
3
4
b, a = sg.butter(1, [0.038, 0.042], 'bandpass') # 9.5Hz - 10.5Hz
sns = sg.filtfilt(b, a, snq)
snq = sns*4+snq
draw(x, snq)

  1. PSD(相敏检波器)处理输入与参考信号(本质为乘法器)
1
2
snqx = snq*square
draw(x, snqx)

1
FFT(snqx, 500)

  1. 低通滤波,减少其他信号对直流分量的影响(等效噪声宽度很小)
1
2
3
b, a = sg.butter(3, 0.00004, 'lowpass') # 0 - 0.01Hz
sns = sg.filtfilt(b, a, snqx)
draw(x, sns)

  1. 直流放大×5
1
2
sns = sns*5
draw(x, sns)