import numpy as np
import matplotlib.pyplot as plt
import math
import sounddevice as sd
import soundfile as sf
import librosa
import os
import wave
# load signals
speech, speech_fs=sf.read("speech.wav")
speech=speech[:,0]
music, music_fs=sf.read("music.wav")
ir, ir_fs=sf.read("impulse_response.wav")
# normalize signals
def normalize(sig):
return sig/max(abs(sig))
speech=normalize(speech)
music=normalize(music)
ir=normalize(ir)
def plotsig(sig, fsr, name, num):
p=plt.subplot(3,1,num)
t=np.linspace(0, len(sig) / fsr, len(sig))
p.plot(t,sig)
p.set_title(name)
p.set_xlabel("$t$ (s)")
p.set_ylabel("Amplitude")
p.grid("both")
plt.figure(figsize=(7,10))
plotsig(speech, speech_fs, "Speech",1)
plotsig(music, music_fs, "Music",2)
plotsig(ir, ir_fs, "Impulse Response",3)
plt.tight_layout()
plt.show()
# Generate signals
fs=16000
t=np.linspace(0,1,fs)
one_khz=0.5*np.sin(1000*2*np.pi*t)
four_khz=0.5*np.sin(4000*2*np.pi*t)
fifteen_khz=0.5*np.sin(15000*2*np.pi*t)
def plotsig(sig, name, num):
p=plt.subplot(3,1,num)
t=np.linspace(0, len(sig) / fs, len(sig))
p.plot(t,sig)
p.set_title(name)
p.set_xlabel("$t$ (s)")
p.set_ylabel("Amplitude")
p.grid("both")
# Plot first millisecond
plt.figure(figsize=(7,10))
plotsig(one_khz[:16],"1 kHz",1)
plotsig(four_khz[:16],"4 kHz",2)
plotsig(fifteen_khz[:16],"15 kHz",3)
plt.tight_layout()
plt.show()
# Play signals (remove comment at the desired signal)
#sd.play(one_khz,fs)
#sd.play(four_khz,fs)
#sd.play(fifteen_khz,fs)
# Generate Signals
fifteen_khz_fs16=fifteen_khz # has already been generated
fs192=192000
t2=np.linspace(0,1,fs192)
fifteen_khz_fs192=0.5*np.sin(15000*2*np.pi*t2)
plt.figure()
plt.plot(t[:16],fifteen_khz_fs16[:16],label="15kHz, fs=16kHz")
plt.plot(t2[:192],fifteen_khz_fs192[:192],label="15kHz, fs=192kHz")
plt.legend()
plt.grid()
plt.xlabel("$t$ (s)")
plt.ylabel("Amplitude")
#for causal system y[n] the impulse response h[n] looks as follows
h = [.2,.2,.2,.2,.2]
def convolve(a,b):
out=np.zeros(len(a)+len(b)-1)
g=np.append(b,np.zeros(len(a)))
x=(a)
for n in range(len(out)):
sum=0
for k in range(len(x)):
sum+=x[k]*g[n-k]
out[n]=sum
return out
a=[1,-1]
b=[2,0,0,1]
print(convolve(a,b))
[ 2. -2. 0. 1. -1.]
# Signals present from above, just need to be shortened to 16 samples
x1=one_khz[:16]
x2=four_khz[:16]
one=one_khz[:16]
four=four_khz[:16]
x1c=convolve(x1,h)
x2c=convolve(x2,h)
plt.figure(figsize=(7,6.6))
plotsig(x1c,"x1 gefaltet mit IR aus 5a",1)
plotsig(one,"x1 gefaltet mit IR aus 5a",1)
plotsig(x2c,"x2 gefaltet mit IR aus 5a",2)
plotsig(four,"x2 gefaltet mit IR aus 5a",2)
plt.tight_layout()
plt.show()
# x1[n]h[n] = h[n]x1[n]
a = convolve(x1,h)
b = convolve(h,x1)
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(10,3))
ax1.stem(a)
ax1.set_title("x1[n]*h[n]")
ax1.set_xlabel("$n$")
ax1.set_ylabel("$amplitude$")
ax2.stem(b)
ax2.set_title("h[n]*x1[n]")
ax2.set_xlabel("$n$")
ax2.set_ylabel("$amplitude$")
plt.tight_layout()
#fig.savefig('homework01_commutativity01.png',dpi=150)
if a.all() == b.all():
print("x1[n]*h[n] equals h[n]*x1[n]")
else:
print("x1[n]*h[n] does not equal h[n]*x1[n]")
x1[n]*h[n] equals h[n]*x1[n]
# (x1[n] h[n]) + (x2[n] h[n]) = (x1[n] + x2[n]) h[n]
a = convolve(x1,h)
b = convolve(x2,h)
c = a+b
d = convolve((x1+x2),h)
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(10,3))
ax1.stem(c)
ax1.set_title("(x1[n]*h[n]) + (x2[n]*h[n])")
ax1.set_xlabel("$n$")
ax1.set_ylabel("$amplitude$")
ax2.stem(d)
ax2.set_title("(x1[n] + x2[n])*h[n]")
ax2.set_xlabel("$n$")
ax2.set_ylabel("$amplitude$")
plt.tight_layout()
#fig.savefig('homework01_commutativity02.png',dpi=150)
if c.all() == d.all():
print("(x1[n]*h[n]) + (x2[n]*h[n]) equals (x1[n] + x2[n])*h[n]")
else:
print("(x1[n]*h[n]) + (x2[n]*h[n]) does not equal (x1[n] + x2[n])*h[n]")
(x1[n]*h[n]) + (x2[n]*h[n]) equals (x1[n] + x2[n])*h[n]