When performing frequency domain (FFT) based processing it is often useful to display a spectrogram of the frequency domain results. While there is a very good SciPy spectrogram function, this takes time domain data and does all of the clever stuff. However if you are processing data in the frequency domain you often just want to build the spectrogram dataset and keep appending FFT results to it.
A spectrogram is a 3D plot, with the following configuration:
- Time is on the X axis
- Frequency is on the Y axis
- Frequency magnitude is shown in colour
This program will use the Matplotlib function imshow() to display the spectrogram.
There are two key tricks to using imshow() for this purpose:
- Rotate the dataset so that time is on the x-axis and frequency is on the y-axis
- Scale the x and y axes labels to correctly show time and frequency
- Remove the second half of the FFT results - once the magnitude of the FFT result has been calculated, the two halves of the result are mirror images so we can discard the upper half.
Here's the code:
import matplotlib.pyplot as pltimport numpy as npfrom scipy import signal# Global ConfigurationplotXAxisSecondsFlag = True # Set to True to Plot time in seconds, False to plot time in samplesplotYAxisHzFlag = True # Set to True to Plot frequency in Hz, False to plot frequency in binsFs = 10000 # Sampling Frequency (Hz)timePeriod = 10 # Time period in secondssampleLength = Fs*timePeriodsinusoidFrequency = 1000 # Frequency of sine wave (Hz)fftLength = 256 # Length of the FFThalfFftLength = fftLength >> 1window = np.hanning(fftLength)time = np.arange(sampleLength) / float(Fs) # Generate sinusoid + harmonic with half the magnitudex = np.sin(2*np.pi*sinusoidFrequency*time) * ((time[-1] - time)/1000) # Decreases in amplitude over timex += 0.5 * np.sin(2*2*np.pi*sinusoidFrequency*time) * (time/1000) # Increases in amplitude over time# Add FFT frames to the spectrogram list - Note, we use a Python list here becasue it is very easy to append tospectrogramDataset = []i = 0while i < (len(x) - fftLength): # Step through whole datasetx_discrete = x[i:i + fftLength] # Extract time domain framex_discrete = x_discrete * window # Apply window functionx_frequency = np.abs(np.fft.fft(x_discrete)) # Perform FFTx_frequency = x_frequency[:halfFftLength] # Remove the redundant second half of the FFT resultspectrogramDataset.append(x_frequency) # Append frequency response to spectrogram dataseti = i + fftLength# Plot the spectrogramspectrogramDataset = np.asarray(spectrogramDataset) # Convert to Numpy array then rotate and flip the datasetspectrogramDataset = np.rot90(spectrogramDataset)z_min = np.min(spectrogramDataset)z_max = np.max(spectrogramDataset)plt.figure()plt.imshow(spectrogramDataset, cmap='gnuplot2', vmin = z_min, vmax = z_max, interpolation='nearest', aspect='auto')plt.title('Spectrogram')freqbins, timebins = np.shape(spectrogramDataset)xlocs = np.float32(np.linspace(0, timebins-1, 8))if plotXAxisSecondsFlag == True:plt.xticks(xlocs, ["%.02f" % (i*spectrogramDataset.shape[1]*fftLength/(timebins*Fs)) for i in xlocs]) # X axis is time (seconds)plt.xlabel('Time (s)')else:plt.xticks(xlocs, ["%.02f" % (i*spectrogramDataset.shape[1]*fftLength/timebins) for i in xlocs]) # X axis is samplesplt.xlabel('Time (Samples)')ylocs = np.int16(np.round(np.linspace(0, freqbins-1, 11)))if (plotYAxisHzFlag == True):plt.yticks(ylocs, ["%.02f" % (((halfFftLength-i-1)*Fs)/fftLength) for i in ylocs]) # Y axis is Hzplt.ylabel('Frequency (Hz)')else:plt.yticks(ylocs, ["%d" % int(halfFftLength-i-1) for i in ylocs]) # Y axis is Binsplt.ylabel('Frequency (Bins)')plt.show()
