Plot DDM, SDM and RF of ILS signal logged with RTL-SDR during ground crossover in front of runway

demod-ils.py

# Plot DDM, SDM and RF of an ILS IF-signal (saved as wav file by SDR#)
# Bandwidth = 32 kHz when input sample rate = 2.048 MHz (decimated by 64)
# Asgaut Eng/2016-10-03

import numpy as np
import scipy.io.wavfile as wav
import scipy.signal as signal
from scipy.fftpack import fft, fftfreq, fftshift
import matplotlib.pyplot as plt

fs, data = wav.read("SDRSharp_20160927_125449Z_111950000Hz_IQ.wav", mmap=True)

T = 1/fs # sample spacing
N = int(0.1 / T) # FFT bin size = 10 Hz

print("Sample rate %f MS/s" % (fs / 1e6))
print("Num samples %d" % len(data))
print("Duration: %f s" % (T * len(data)))
iq = np.empty(N, np.complex64)

Nd = N // 8 // 8
Td = T * 8 * 8
print("Decimated sample rate: %f S/s" % (1 / Td))

t = np.linspace(0, 90, 900) # from 0 to 90 s in 900 steps
ddm = np.empty(len(t), np.float)
sdm = np.empty(len(t), np.float)
rf = np.empty(len(t), np.float)
for i in range(0, len(t), 1):
    index = int(t[i] * fs)
    iq.real, iq.imag = data[index:index+N,0], data[index:index+N,1]
    iqd = signal.decimate(iq, 8, zero_phase=False)
    iqd = signal.decimate(iqd, 8, zero_phase=False)
    # Demodulate AM signal
    demod = np.abs(iqd)
    f = 1/Nd * fft(demod)
    dc = np.abs(f[0])
    # Get 90 and 150 modulation levels (assumes FFT bin size = 10 Hz)
    mod90 = (abs(f[90//10]) + abs(f[len(f)-90//10]))/dc*100
    mod150 = (abs(f[150//10]) + abs(f[len(f)-150//10]))/dc*100
    ddm[i] = mod150 - mod90
    sdm[i] = mod150 + mod90
    rf[i] = dc

plt.subplot(311)
plt.title("DDM")
plt.plot(t, ddm)
plt.grid()
plt.subplot(312)
plt.title("SDM")
plt.plot(t, sdm)
plt.grid()
plt.subplot(313)
plt.title("RF level")
plt.plot(t, rf)
plt.grid()
plt.show()

Thank you.
Can you please provide me the wav file and the measured ddm, sdm, rf level with PIR?