Thursday, September 8, 2016

Homemade AM modulator

I recently got into Python and wrote an AM modulator. I learned the idea from a free PDF called Think DSP. The concept is really simple and just involves multiplying audio samples by sine wave samples. It was a bit difficult where I had to cast bytes into floats and then back again. This is a no-brainer in Liberty BASIC, but I stuck with Python and got it working the same day.

I have previously tried to write processing software for sound but failed. I thought the sound samples might be a weird logarithm that must be undone, but during this project I saw that sound is intuitive to work with and produces the expected results if you use 8-bit unsigned WAV. Being unsigned, it gives you values from 0 to 255, which really is simple to manipulate, and Audacity will show the expected waveform when you open it.

My Python script takes an input wave file and multiplies 100,000 samples by a sine wave. To save the trouble of having to use Audacity's Raw Audio import feature, the script copies the WAV headers to the new file. The output audio length therefore isn't accurate but Audacity and HDSDR have safeguards to avoid reading too far.

Here's the resulting output in HDSDR.
















One more thing before I share the script: this is a mono WAV file, which is why HDSDR's frequency scale starts at 0 on the left. I tested with a stereo 8-bit WAV and HDSDR considers it to be IQ and centers 0 on the frequency scale, with negative and positive frequencies on either side. This is good to know for when I start generating quadrature signals. 


And now the script:
-----------------
#import matplotlib.pyplot as plt
import math

def radians(degrees):
    return (degrees/360)*2*math.pi


#plt.axis([0,1000000,-1,1])
#plt.ylabel('some numbers')
x=[]
y=[]
multiplier=[]
demod=[]
amp=1
freq=18 #18 corresponds to 15.075 kHz in this case (see picture above)
phase=0
for d in range(1000000):
    x.append(d)
    #y.append(amp*math.sin((2*math.pi*freq*radians(d))+phase))
    y.append(math.sin((2*math.pi*freq*radians(d))+phase))
    #multiplier.append(amp*math.sin((2*math.pi*freq*radians(d))+phase))

#for element in range(len(y)):
    #demod.append(y[element]*multiplier[element])
#plt.plot(x,demod)
#plt.plot(x,y)
#plt.show()
bytepos=0
with open("D:/time8a.wav","rb") as i: #Input File
    with open("D:/time_am_hf18.wav","wb") as o: #Output file
        
    #byte = i.read(1)
        for idx in range(0,44):
            byte=i.read(1)
            o.write(byte)
        
        while byte:
            bytepos +=1
            test=float(float(int.from_bytes(byte,byteorder="big"))*((y[bytepos]/2)+0.5))
            #print(int.from_bytes(byte,byteorder="big")," * ",y[bytepos]," = ",float(float(int.from_bytes(byte,byteorder="big"))*y[bytepos]))
            tmp = [int(test),]
            #print(bytes(tmp)," (",int(test),")")
            o.write(bytes(tmp))
            byte = i.read(1)
-----------------

No comments:

Post a Comment