python3で『FM変調』と『ダイレクトコンバージョン』

WFM変復調。IQ復調とFIRフィルターを使用。

import numpy as np

from scipy import signal as sg

from scipy.fftpack import fft,ifft,fftshift,fftfreq

import matplotlib.pyplot as plt

fb = 440 #baseband

fc = 10000 #carrier

fs = 44100 #sampling freq

A = 1 #amp

t = 3 #time

m = 10 #Frequency shift m > 1 wideFM

dat = np.arange(0,t,1/fs)

pt = 2*np.pi*dat

#sin(alpha)cos(beta)

sig = np.sin(fb * pt)

carrier = np.cos(fc * pt)

#

wave = A * np.cos((fc*pt) + m * np.sin(fb*pt))

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

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

#plt.plot(wave[0:500])

#plt.show()

#解析的信号 hilbert transform 

'''

#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(wave)

asig = wave + (hil * 1j)

Q=np.real(asig)

I=np.imag(asig)

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

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

#plt.show()

amp = np.abs(asig)

# under the same mean

phase = np.unwrap(np.angle(asig))

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

#demodulation

want= np.sin(fb * pt)

# normalize

nyq = fs/2 

Nf = (freq - fc)/ nyq 

#low pass filter

# num must be odd number

num = 15

cutoff = (fc/2) / nyq

lpf = sg.firwin(num,cutoff)

signal = sg.lfilter(lpf,1,Nf)

#　Amplitude correction

signal = (signal * 5)

plt.plot(signal[0:500])

plt.plot( want[0:500])

plt.show()