rtl_fm_analysis.py

#!/usr/bin/env python

from pylab import *
from scipy import signal
import struct
import sys
import numpy

def plotSpectrum(y, Fs, titleStr):
    """
    Plots a Single-Sided Amplitude Spectrum of y(t)
    """
    figure()
    n = len(y) # length of the signal
    k = arange(n)
    T = n/Fs
    frq = k/T # two sides frequency range
    frq = frq[range(n/2)] # one side frequency range

    Y = numpy.fft.fft(y)/n # fft computing and normalization
    Y = Y[range(n/2)]
    print "*"*10
    print "The aveage "
    avgFs = (sum(abs(Y)))/float(len(Y))
    print avgFs
    print "*"*10

    plot(frq,abs(Y),'r') # plotting the spectrum
    title(titleStr)
    xlabel('Freq (Hz)')
    ylabel('|Y(freq)|')


N = 256
Wc = 0.208
order = 3
f = open("audio.raw", "rb")

Fs = 48000.0
i = 0
dataSet = []
while i != 256:
    #try:
    dat = f.read(2)
    (temp, )= struct.unpack('h', dat)
    dataSet.append(temp)
    i = i+1
    #finally:
    #    print i
    #    print "Out of inputs"
    #    f.close

plotSpectrum(dataSet, Fs, 'audio')
[b, a] = signal.butter(order, Wc, 'high')

calYn = signal.lfilter(b, a, dataSet)

# Set initial conditions
Yn = [0]*(4+len(dataSet))
Xn = [0, 0, 0, 0] +  dataSet[0:]
# Filter by hand
for i in range(3, len(dataSet)+3):
    Yn[i] = 0.51378993*Xn [i] - 1.54136978*Xn [i-1] + 1.54136978*Xn [i-2] - 0.51378993*Xn [i-3] \
        + 1.71160576*Yn[i-1] - 1.13513329*Yn[i-2] + 0.26358035*Yn[i-3]

Yn = Yn[4:]

print Yn
print calYn

plotSpectrum(calYn, Fs, 'Calculated using lfilter')
plotSpectrum(Yn, Fs, 'Calculated by implementing the DE')
show()