I have a recording device that exhibits Wow and Flutter, or periodic frequency deviation. I use this device to record a calibrated 3.15kHz test tone from a tape film to generate a .wav file on my PC.

In an ideal world, an FFT on the data contained in the .wav file would result in a single sharp peak at 3.15kHz. In the real world, there is deviation in this frequency due to mechanical errors such as head scrape, non-concentric film rollers, motor vibration which all cause the tape film to either speed up or slow down slightly.

These deviations are called Wow and Flutter, and appear as errors frequency modulated by the 3.15kHz test tone.

What I would like to do:

Assume I have the .wav file, sampled at 48kHz. Let's say it contains an FM waveform with carrier frequency at 3.15kHz, and a message signal with bandwidth ~200Hz.

I'd like to build a python script to analyse and demodulate the signal to give me just the 200Hz signal. I've tried several models across the internet simulating PLLs, slope detectors etc. But these are all hard coded difficult to understand and I just can't seem to mould them to be used in my application.

I want to try implementing an IQ FM demodulator but I am unable to find a good reference to how FM waves are modulated via IQ signals and I definitely can't find one for demodulation.

Can someone point me in the right direction and let me know off the bat if there are any problems or roadblocks with what I'm trying to do?

  • $\begingroup$ In an ideal world, an FFT on the data contained in the .wav file would result in a single sharp peak at 3.15kHz: Only if your DFT length and sampling rate happen to divide such that there is a bin at exactly that frequency. WLOG, there isn't, and with limited clock accuracy and sample depth and computational precision in the real world, there actually can't be. $\endgroup$ – Marcus Müller Mar 23 '17 at 22:38
  • $\begingroup$ You can't extract the input 200 Hz signal if the frequency variations due to device imperfections happen to fall into the same bandwidth. There's simply mathematically no difference between FM modulation and changing the frequency of a tone. $\endgroup$ – Marcus Müller Mar 23 '17 at 22:39
  • $\begingroup$ @MarcusMüller The 200Hz signal is the frequency variations due to device imperfection. I'm trying to measure the frequency deviation and associated energy levels. $\endgroup$ – Josef de Joanelli Mar 23 '17 at 22:51
  • $\begingroup$ Do our answers that we gave here help you? electronics.stackexchange.com/questions/293706/… $\endgroup$ – Dan Boschen Mar 23 '17 at 23:35
  • $\begingroup$ @DanBoschen Those answers are sort of helpful, but I'm confused about the frequency discriminator, which I'm guessing only pertains to discrete data because it is a keying modulation. My 200Hz signal is a continuous signal that I'm trying to extract. $\endgroup$ – Josef de Joanelli Mar 24 '17 at 0:05

Frequency is the derivative (or 1st difference) of phase. To get an IQ signal from which you can get phase angles, first multiply your wave file samples by a cosine at 3.15kHz and a sinewave at a frequency of 3.15kHz. That will hetrodyne your signal down to baseband IQ. Use atan2 of the IQ array to get an array of angles, take the 1st difference between successive angles, and note that big jumps in phase differences will likely need to be unwrapped to a smaller range. Then low pass filter (with a 200 Hz cutoff) the result. That's a simplified FM demodulator . See if that looks close to what you are looking for.


There are several FM demodulator solutions at this post that many of us here have responded to that may help you: https://electronics.stackexchange.com/questions/293706/fsk-demodulation-using-dsp/293723?noredirect=1#comment675738_293723.

In addition to what is contained in that link, given what you described in that you have I and Q signals, and if you are indeed post-processing data and not trying to make an efficient implementation, I would likely do the following:

Do a single sided down-conversion to remove your carrier by mutliplying by $e^{-j2 \pi 3150 t}$. This will move your modulated signal now centered on 3.15 KHz to DC.

You will want to have an idea of the modulated bandwidth of your signal and then low pass this baseband signal prior to doing the actual FM demodulation (this step is important since the demod process is a high pass that will amplify higher frequency content!). Your actual bandwidth will depend on the modulation index (which is the ratio of the rate of your modulation signal versus the amount the magnitude of the signal causes the frequency to deviate: $\frac{f_{mod}}{f_{dev}}$). The 200 Hz you specified I believe is the rate of your modulated signal, so the question regarding bandwidth is how far in frequency does modulated signal pull your carrier frequency? See Carson's rule for more details on the rough estimate, which is basically adding the rate of modulation to the frequency deviation to get your total bandwidth. Once you establish your bandwidth either from this approach or inspection of the signal spectrum, you want to low pass filter the down-converted signal to reasonably close to this bandwidth (single-sided) to minimize out of band distortion from effecting the measurement.

Then check that you don't have a significant residual frequency offset (because your carrier may not have been exactly as where you thought it was) by looking at the unwrapped average phase of your signal, using Phase = ATAN2(Q,I), and if in Matlab unwrap(phase). Take a moving average so that you can observe the trend- Your frequency error (since we could not down-convert with your exact carrier unless you were synchronized) will be the slope of this trend line ($\frac{d\phi}{dt}$). Remove the bulk of any residual frequency error using the same exponential rotation we used to down-convert (just in this case it will be a much smaller frequency offset and you will need to go either positive or negative depending on the direction of the trendline. The goal is to visibly flatten the general trend of the phase but no need to be entirely flat as any residual slope will just be a frequency offset, which means a DC offset in your final demodulated signal.

Once you have the trend line reasonably flattened, you can convert this phase versus time to frequency versus time by simply taking a derivative (digitally difference successive samples), and then low pass your final result.

  • $\begingroup$ What do you mean unwrap the phase? $\endgroup$ – Josef de Joanelli Mar 24 '17 at 5:39
  • $\begingroup$ Are you using Matlab? It has a command unwrap(phase) as otherwise the phase will roll at the $2\pi$ boundaries. If the phase had an upward slope (for example) if it was plotted normally you would see a sawtooth as it would not go past $2\pi$. By unwrapping it you end up getting a single smooth function of phase growing continuously with time. Since $d\phi/dt$ is frequency, the slope of that line would be frequency (when divided by $2\pi$), so that is what you would want to work with. Make sense? $\endgroup$ – Dan Boschen Mar 24 '17 at 5:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.