ECE3340 Numerical Methods

Jupyter notebook (translated from Mathematica)

In support of HW 4A¶

ECE3340 - University of Houston

Han Q. Le copyrighted¶

In the follow, we will try to follow as closely to the Mathematica code as possible (see this and this), in other words, "literal translation," even if it might not be optimal with Python, because we want to help those getting start with Python and Jupyter notebook to have the same sense of familiarity as with Mathematica notebook, and vice versa.

In [1]:
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import sounddevice as sd

1. A binary digital message¶

In [2]:
msg = [1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0
   , 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1]
In [3]:
# while a simpler function can be written, this is trying to do literal translation
# from Mathematica for illustration
def  unitbox(t,m):
    if np.abs(t-m)< 0.25 :
        r=1
    else:
        r=0
    return r
def digitsig(t,messg):
    s=sum([messg[m-1]*unitbox(t,m) for m in np.arange(len(messg))+1])
    return s
In [4]:
srate=8 ; deltt=1./srate
tarr=np.arange(0,len(msg)+1,deltt)
sig=[digitsig(u,msg) for u in tarr]
sigplot=plt.plot(sig)
plt.show(sigplot)
In [5]:
deltF=1./(len(tarr)*deltt); npt=len(sig)
farr=(np.arange(npt)-1-np.floor(npt/2))*deltF
plt.plot(farr)
Out[5]:
[<matplotlib.lines.Line2D at 0x1e0a81787b8>]

We will stick to Mathematica convention of Fourier (physics based), hence we use ifft (inverse fft) of numpy so that it is consistent with our Mathematica Tutorial results. In addition, we include the constant factor.

In [6]:
sigFT=np.fft.ifft(sig)*(deltt/np.sqrt(2*np.pi))
sigFTsh=np.fft.fftshift(sigFT)
sigFTAbs=np.absolute(sigFTsh)
sigFTdB=20*np.log10(sigFTAbs)

C:\ProgramData\Anaconda3\lib\site-packages\ipykernel_launcher.py:4: RuntimeWarning: divide by zero encountered in log10
  after removing the cwd from sys.path.
In [7]:
plt.plot(farr,sigFTdB)
Out[7]:
[<matplotlib.lines.Line2D at 0x1e0a87b45c0>]

2. Signal with carrier¶

In [4]:
def carrier(fc,t):
    c=np.cos(2*np.pi*fc*t)
    return c
In [5]:
srate2 = 1100; deltt2 = 1./srate2;
tarr2 = np.arange(0,len(msg)+1,deltt2)
sig2 = np.array([digitsig(u,msg) for u in tarr2])
fc = 125.
carr2 =np.array(carrier(fc, tarr2))
sigcarr = sig2*carr2
imgsz= plt.figure(figsize=(8,8*0.1));sigcarrplt=plt.plot(sigcarr)
plt.show(sigcarrplt)
In [10]:
sndsr=22050;
sd.play(sigcarr,sndsr)
In [6]:
deltF=1./(len(tarr2)*deltt2); npt2=len(sigcarr)
farr2=(np.arange(npt2)-1-np.floor(npt2/2))*deltF
plt.plot(farr2)
Out[6]:
[<matplotlib.lines.Line2D at 0x15027f330b8>]
In [12]:
sigcarrFT=np.fft.ifft(sigcarr)*(deltt2/np.sqrt(2*np.pi))
sigcarrFTsh=np.fft.fftshift(sigcarrFT)
sigcarrFTAbs=np.absolute(sigcarrFTsh)
sigcarrFTdB=20*np.log10(sigcarrFTAbs)
specSC=plt.plot(farr2,sigcarrFTdB)
plt.show(specSC)
In [13]:
sigcarrdBmax=np.ma.max(sigcarrFTdB)
plt.xlim(fc-5,fc+5);plt.ylim(sigcarrdBmax-60,sigcarrdBmax+5)
plt.plot(farr2,sigcarrFTdB)
Out[13]:
[<matplotlib.lines.Line2D at 0x1e0a9034860>]

1.3 Signal with carrier + noise¶

In [7]:
noise = np.random.normal(0.,0.5,len(sigcarr));
# this is very bad noise by normal comm stardard
sigcarrnoise = sigcarr + noise;
imgsz= plt.figure(figsize=(8,8*0.1));plt.plot(sigcarrnoise)
Out[7]:
[<matplotlib.lines.Line2D at 0x15027edb5c0>]
In [15]:
sndsr=22050;
sd.play(sigcarrnoise,sndsr)
In [8]:
sigcarrNFT=np.fft.ifft(sigcarrnoise)*(deltt2/np.sqrt(2*np.pi))
sigcarrNFTsh=np.fft.fftshift(sigcarrNFT)
sigcarrNFTAbs=np.absolute(sigcarrNFTsh)
sigcarrNFTdB=20*np.log10(sigcarrNFTAbs)
sigcarrNdBmax=np.ma.max(sigcarrNFTdB)
plt.xlim(fc-20,fc+20);plt.ylim(sigcarrNdBmax-60,sigcarrNdBmax+5)
specSCN=plt.plot(farr2,sigcarrNFTdB)
plt.show(specSCN)

1.4 Band-pass filter¶

In [9]:
sigoutArr=[]; sigFTfiltArr=[]
for bw in [5,10,20,30]:
    filtBP=[1 if np.absolute(np.absolute(f)-fc) <= 0.5*bw else 0 for f in farr2]
    sigFTfiltered0 = sigcarrNFTsh*np.array(filtBP);
#    plt.xlim(fc-20,fc+20)
    sigFTfiltArr.append(sigFTfiltered0)
    sigFTfilt1 = np.fft.ifftshift(sigFTfiltered0);
    sigout=np.real(np.fft.fft(sigFTfilt1))*deltF*np.sqrt(2*np.pi)*len(sigFTfilt1)
    plt.figure(figsize=(4,3)); plt.xlim(0,3)
    plt.plot(tarr2,sigout)
    sigoutArr.append(sigout)

Observe that we include those constant factors for fft and ifft only to show that the signal amplitude is back where it should be. In practice, where amplitude doesn't matter, the factors are often neglected.

In [17]:
for s in sigFTfiltArr:
    plt.figure(figsize=(6,3)); plt.xlim(fc-20,fc+20); 
    plt.plot(farr2,10*np.log10(np.absolute(s)))

C:\ProgramData\Anaconda3\lib\site-packages\ipykernel_launcher.py:3: RuntimeWarning: divide by zero encountered in log10
  This is separate from the ipykernel package so we can avoid doing imports until

1.5 Remove the carrier¶

although we already had the filtered result with bw=30 in sigFTfiltArr list, we can do it again for reinforced learning

In [18]:
bwbp = 30;
filtBP=[1 if np.absolute(np.absolute(f)-fc) <= 0.5*bw else 0 for f in farr2]
sigFTfiltered0 = sigcarrNFTsh*np.array(filtBP);
sigFTfilt1 = np.fft.ifftshift(sigFTfiltered0);
sigout=np.real(np.fft.fft(sigFTfilt1))*deltF*np.sqrt(2*np.pi)*len(sigFTfilt1)
plt.figure(figsize=(8,0.15*8))
plt.plot(sigout)
Out[18]:
[<matplotlib.lines.Line2D at 0x1e0aa35cac8>]

the above should be compared with the noisy signal before filtering

image.png

In [18]:
sigReceived=[]; 
sigP=sigout**2; # power
sigoutPFT=np.fft.fft(sigP)
sigoutPFTsh=np.fft.fftshift(sigoutPFT)*deltt2/np.sqrt(2*np.pi)
for bwlp in [2,5,10,20] :
    filtLP=[1 if np.absolute(f) <= bwlp else 0 for f in farr2]
    sigFTfiltLP= np.fft.ifftshift(sigoutPFTsh*np.array(filtLP));
    sigoutBB=np.real(np.fft.ifft(sigFTfiltLP))*deltF*np.sqrt(2*np.pi)*len(sigFTfiltLP)
    plt.figure(figsize=(8,0.15*8))
    sigReceived.append(plt.plot(sigoutBB))

How to use Fourier array shift functions:
fftshift -> RotateRight
ifftshift -> RotateLeft
examples below

In [2]:
a =np.fft.fftshift(range(11))
print(a)

[ 6  7  8  9 10  0  1  2  3  4  5]
In [5]:
np.fft.ifftshift(a)
Out[5]:
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10])
In [ ]: