日常雑記2020年7月28日＿Python_AM_modem

暇なので『AM変調』と『ダイレクトコンバージョン』による復調

処理をPythonで処理中。以下は『AM変調』後に『ノイズ』を付加

したものを復調。環境はLinuxでPython3。誰かの参考になれば。

import numpy as np
from scipy import signal as sg
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt

fb = 440 #baseband
fc = 10000 #carrier
fs = 44100 #sampling freq
A = 1 #amp
t = 1 #time
m = 0.5 #Modulation index

dat = np.arange(0,t,1/fs)
pt = 2*np.pi*dat

#sin(alpha)cos(beta)
sig = A + m * np.sin(fb*pt)
carrier = A * np.cos(fc*pt)

#under the same mean
wave = sig * carrier
#wave = carrier + ( A*np.sin(pt*(fb+fc)) + A*np.sin(pt*(fb-fc)) )/2 * m

#plt.plot(sig[0:100])
#plt.plot(carrier[0:100])
#plt.plot(wave[0:100])
#plt.show()

#解析的信号 hilbert transform
'''
# hilbert transform4 steps
Num=len(wave)
X=fft(wave,Num)
Z=np.hstack((X[0],2*X[1:Num//2],X[Num//2],np.zeros(Num//2-1)))
asig=ifft(Z,Num)
'''

#AWGN noise add
#Add AWGN noise to the transmitted signal
nMean = 0 #noise mean
nSigma = 0.1 #noise sigma
n = np.random.normal(nMean, nSigma, len(dat))
noisy = wave + n  #noisy received signal

#plt.plot(wave[0:100])
#plt.plot(noisy[0:100])
#plt.show()

#analytic signal
hil = sg.hilbert(noisy)
asig = noisy + (hil * 1j)
Q=np.real(asig)
I=np.imag(asig)

#plt.plot(I[0:100])
#plt.plot(Q[0:100])
#plt.show()

# under the same mean amplitude
#envelope = np.sqrt(I*I+Q*Q)
envelope = np.abs(asig)

#plt.plot(envelope[0:100])
#plt.show()

amp = envelope
# under the same mean
#phase=np.arctan2(np.imag(asig),np.real(asig))*180/np.pi
phase = np.unwrap(np.angle(asig))

#freq = np.diff(phase)/(2*np.pi)*fs

# A*e(phi *1j) signal
aej=(amp) * np.exp(phase * 1j)

#変調処理の逆を行う
dband = aej - carrier
want=np.abs(sg.hilbert(np.imag(dband)))
#欲しい信号
plt.plot(want[0:300])
plt.plot(sig[0:300])
plt.show()