Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
PSK63 IMD measurement
#!/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')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment