仿真报告1:模拟通信系统抗噪声性能仿真

调制解调器基本原理图

调制器

调制器信号时域分析

令$$m(t)=cos(2\pi f_{m}t),f_m=10Hz,$$
$$c(t)=cos(2\pi f_{c}t),f_c=100Hz$$
则通过调制器得到$s(t)$信号为
$$s(t)=cos(2\pi f_{m}t)cos(2\pi f_{c}t)$$

from matplotlib import pyplot as plt
import numpy as np
from scipy.fftpack import fft,ifft
from matplotlib.pylab import mpl

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

#定义m(t)、s(t)
t = np.linspace(-np.pi/6,np.pi/6,1024,endpoint=True)
fm = 10 #Hz
mt = np.cos(2*np.pi*fm*t)
fc = 100 #Hz
ct = np.cos(2*np.pi*fc*t)
plt.subplot(421)#画c(t)
plt.plot(t,ct)
plt.xlabel("Time t (s)")
plt.ylabel("振幅 c(t)")
plt.subplot(423)
plt.plot(t,mt)
plt.xlabel("Time t (s)")
plt.ylabel("振幅 m(t)")
#plt.show()

运行得到$m(t)、c(t)、s(t)$时域波形

快速傅里叶变换FFT

已知了时域信号
$$m(t)=cos(2\pi f_{m}t),f_m=10Hz,$$
$$c(t)=cos(2\pi f_{c}t),f_c=100Hz$$
其相应频域信号
$$m(f)=\frac{1}{2}[\delta (f-f_m)+\delta (f+f_m)]$$
$$c(f)=\frac{1}{2}[\delta (f-f_c)+\delta (f+f_c)]$$
由于$s(t)$为$m(t)、c(t)$的时域乘积，则$s(f)$为$m(f)、c(f)$的频域卷积。通过冲激函数的频移性质易得
$$\begin{gathered} s(f)=\frac{1}{4}[\delta (f-f_m-f_c)+\delta (f-f_m+f_c) \\ +\delta (f+f_m-f_c)+\delta (f+f_m+f_c)] \end{gathered}$$
程序验证

#c(f)
cf = fft(ct)
N = len(cf)
f = np.arange(N)
abs_cf = np.abs(cf)
abs_cf = abs_cf/N #归一化
f = f[range(int(N/2))]#单边化
abs_cf = abs_cf[range(int(N/2))]#单边化
plt.subplot(422)
plt.plot(f,abs_cf)
plt.xlabel("Frequency f (hz)")
plt.ylabel("振幅 c(f)")

#m(f)
mf = fft(mt)
N = len(mf)
f = np.arange(N)
abs_mf = np.abs(mf)
abs_mf = abs_mf/N #归一化
f = f[range(int(N/2))]#单边化
abs_mf = abs_mf[range(int(N/2))]#单边化
plt.subplot(424)
plt.plot(f,abs_mf)
plt.xlabel("Frequency f (hz)")
plt.ylabel("振幅 m(f)")

#s(t) = m(t)×s(t)
st = mt*ct
plt.subplot(425)
plt.plot(t,st)
plt.xlabel("Time t (s)")
plt.ylabel("振幅 s(t)")
#plt.show()

#对s(t)进行傅里叶变换，s(f) = m(f)*s(f)
sf = fft(st)
N = len(sf)
f = np.arange(N)
abs_sf = np.abs(sf)
abs_sf = abs_sf/N #归一化
#单边化
f = f[range(int(N/2))]
abs_sf = abs_sf[range(int(N/2))]
plt.subplot(426)
plt.plot(f,abs_sf)
plt.xlabel("Frequency f (hz)")
plt.ylabel("振幅 s(f)")
#plt.show()

频域波形

解调器

解调部分

$$s_d(t)=s(t)\cdot c(t)$$

sdt = st*ct
plt.subplot(427)
plt.plot(t,sdt)
plt.xlabel("Time t (s)")
plt.ylabel("振幅 Sd(t)")

$$s_d(t)\stackrel{FFT}{⟶}{s}_d(f)$$

#Sd(f)
sdf = fft(sdt)
N = len(sdf)
f = np.arange(N)
abs_sdf = np.abs(sdf)
abs_sdf = abs_sdf/N
#单边化
f = f[range(int(N/2))]
abs_sdf = abs_sdf[range(int(N/2))]
plt.subplot(428)
plt.plot(f,abs_sdf)
plt.xlabel("Frequency f (Hz)")
plt.ylabel("振幅 Sd(f)")

plt.tight_layout()
#plt.show()

$s_d(t)、s_d(f)$波形

滤波部分——h(t)设计

还记得我们最开始定义的$m(t)$吗？$m(t)$的频率$f_m$为10Hz
通过这一步我们要把$s_d(f)$中属于$m(t)$的那一部分频率给滤出来
观察$s_d(f)$波形，可以发现我们其实只需要把波形中10Hz出的峰峰给保留下来，其他的峰峰都给置零就行了

从而我们只需要用一个合适的矩形信号$h(f)$在频域和$s_d(f)$相乘就可以了

频域矩形信号对应时域抽样信号
$$Rect(f)\stackrel{ifft}{⟶}Sa(t)=\frac{sin(wt)}{wt}$$

#h(t)
ht = 30*(np.sin(30*np.pi*t)/(30*np.pi*t))
plt.figure()
plt.subplot(221)
plt.plot(t,ht)
plt.xlabel("Time t (s)")
plt.ylabel("振幅 h(t)")

#h(f)
hf = fft(ht)
N = len(hf)
f = np.arange(N)
abs_hf = np.abs(hf)
abs_hf = abs_hf/N #归一化
#单边化
f = f[range(int(N/2))]
abs_hf = abs_hf[range(int(N/2))]
plt.subplot(222)
plt.plot(f,abs_hf)
plt.xlabel("Frequency f (Hz)")
plt.ylabel("振幅 h(f)")

#解调：mo(f) = sd(f)×h(f)
mof = sdf * hf
f = np.arange(len(mof))
plt.subplot(223)
f1 = f[range(int(N/2))]#单边化
abs_mof = np.abs(mof)[range(int(N/2))]/len(f)#单边归一化
plt.plot(f1,abs_mof)
plt.xlabel("Frequency f (Hz)")
plt.ylabel("振幅 mo(f)")

#mot
mot = ifft(mof)/len(f)*2
N = len(mot)
t = np.arange(N)/len(f)
t = t[range(int(N/2))]#单边化
mot = mot[range(int(N/2))]#单边化
plt.subplot(224)
plt.plot(t,mot)
plt.xlabel("Time t (s)")
plt.ylabel("振幅 mo(t)")

plt.tight_layout()
plt.show()

结论与分析

原始信号$m(t)$与调制解调信号$m_o(t)$比对

$m(t)$
$m_o(t)$

上图可见$m(t)$与$m_o(t)$时域波形一致，继续增加时域抽样点的个数可以进一步提高解调信号的还原精度

频域信号$m(f)$与$m_o(f)$，可z见二者频域波形一致，均为10Hz