PSK63 IMD measurement

imd.py

#!/usr/bin/env python3

import scipy.io.wavfile
from scipy.fftpack import fft, fftfreq
from scipy.signal import flattop
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import numpy as np
import sys

sample_rate, samples = scipy.io.wavfile.read(sys.argv[1])

N = 4096

w = flattop(N)
norm_w = np.sqrt(np.sum(np.square(w))/N)

centre_freq = 1500
width = 1000
bin_width = sample_rate/N
b_m = int((centre_freq - width/2)/bin_width)
b_p = int((centre_freq + width/2)/bin_width)
xf = fftfreq(N, 1.0/sample_rate)
yf = fft(samples[:N] * (w/norm_w) / np.iinfo(np.int16).max)
yfdB = 20*np.log10(np.abs(yf[b_m:b_p])/N) + 10*np.log10(2.0)

callsign, freq, time = sys.argv[1].split('-')

fig, ax = plt.subplots(1,1)
ax.plot(xf[b_m:b_p], yfdB)
ax.yaxis.set_major_locator(ticker.MultipleLocator(5))
plt.grid()
plt.xlabel('Audio Frequency (Hz)')
plt.ylabel('dBFS')
plt.title('{}    {:.3f}kHz    {}:{}:{}'.format(callsign.upper(), int(freq)*1e-3, time[:2], time[2:4], time[4:6]))
#plt.show()
plt.savefig(sys.argv[1][:-3] + 'png')