学习笔记 | Python模拟锁相放大器的工作过程

发表于2021-05-07更新于2022-06-23字数统计1.5k阅读次数

准备

本篇笔记将使用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=0∞sin[(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=C∫A2dtW=\int Pdt=C\int A^2dt代替，因为信号和噪声的积分时间是相同的。而信号的波形是采样得到的，所以WW可以用∑xi2\sum x_i^2代替，于是有：

SNR=10lg∑si2∑ni2⟹∑ni2=∑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)