3
$\begingroup$

The script essentially pulls PCM data from the sound card and stores them in a buffer. Then I apply fft and graph it. I am looking through this blog post for some examples and ideas on how to implement it myself.

Here is the gist of the script:

self.rate = 48100
self.buffer_size = 1024
...
def get_fft(self):
        data = self.get_pcm() #pulls 1 buffer (1024) of PCM data from the soundcard.
        fft = numpy.abs(numpy.fft.fft(data))
        fft = numpy.divide(fft, 10000)
        freq = numpy.fft.fftfreq(len(fft), 1.0 / self.rate)
        return freq[:int(len(freq)/2)], fft[:int(len(fft)/2)]

I am graphing it in a window with PyQtGraph and using a simple tone generator website to produce frequencies. enter image description here

The frequencies do appear in the correct spot but their magnitude varies quite a bit even though I am not changing the volume or anything. This is especially the case when I slide toward the lower frequencies < 100hZ and the data is much higher in magnitude compared with the rest of the graph.

Am I doing something wrong, and if not, is there a way to dampen the lower frequencies so the entire spectrum appears "equal in magnitude" so to say?

UPDATE 1 (added a window function to the PCM data before FFT:

enter image description here

$\endgroup$
  • $\begingroup$ My first suspicion is that you might be see scalloping loss from locations that are in between the midpoints of the fft bins. @RichardLyons made this nice white paper describing the effect: dspguru.com/sites/dspguru/files/…. If you did a much finer scan over a narrower frequency range to cover a few bins (looks like you have 1024 frequency bins), that could be confirmed or the pattern may suggest other theories. $\endgroup$ – Dan Boschen Feb 26 '17 at 11:28
  • $\begingroup$ @DanBoschen Thanks for the response, this does seem like the issue I'm having. Increasing the buffer size slows down the processing aspect (I'm trying to go for a real-time approach) so I will look at interpolating it as mentioned in the answer below. $\endgroup$ – skyguy126 Feb 26 '17 at 18:13
  • $\begingroup$ Did you try applying a window function first? fft = numpy.abs(numpy.fft.fft(data * numpy.blackman(self.buffer_size))) $\endgroup$ – endolith Feb 27 '17 at 21:43
3
$\begingroup$

A very simple solution to your problem is to use a window on your dataset prior to taking the FFT. This classic paper by fred harris details the scalloping loss for all the popular windows. Note the following from this paper:

The scalloping loss for a rectangular window (meaning no window, which I believe is what you are doing) is 3.92 dB. The minimum listed in Table 1 in his paper is 0.83 dB using a "Minimum 4-sample Blackman Harris". This presents a very simple solution in that you multiply your time domain by the appropriate window function prior to taking the fft, reducing the amplitude variation from 3.92 dB (36.4% drop from peak at bin centers) down to 0.83 dB (9.12% drop)!

Further you can combine this with overlap-add techniques if you are block processing streaming data, which is also further explained in the paper.

https://www.utdallas.edu/~cpb021000/EE%204361/Great%20DSP%20Papers/Harris%20on%20Windows.pdf

Also if your requirements is simply for graphical reasons for a swept sinusoidal signal, do know that the signal is not actually "lost" but distributed to adjacent bins. You can recover most of the signal by root sum-squaring the adjacent bins (or all the signal by rss'ing all the bins and displaying this at the max signal location). Doing this, for a sinusoidal signal, will eliminate all inter-bin variation.

Alternatively similar to the previous paragraph for a streaming application where we want to compensate the output for "presentation" purposes on a graph, predistortion (equalization) of the passband ripple prior to the fft would be simple to do, as long as you have more data available than the size of your FFT. A 3 tap FIR compensator would bring the ripple from 3.9 dB to less than 1 dB. This technique could also be combined with the windowing mentioned to bring the ripple down well below the 0.83 dB provided by the window. Without windowing, each bin of the FFT output is a Sinc filter at the bin's center frequency with the nulls at all the other bins. Recognizing this, we could use an inverse Sinc approach as a compensator, following the design for a 3 tap inverse Sinc compensator given as my response to this post:

how to make CIC compensation filter

enter image description here

enter image description here

Trying out a first pass of this approach on a simple 4 pt DFT (to match the plots above), I was able to achieve the result shown in the plot below showing the scalloping loss of -3.92 dB as listed in fred harris' paper, and then the reduced deviation after using a simple 3 tap FIR prefilter. This can be reduced further with an IIR approach of similar order, or more taps but at that point I would likely use a windowing approach (as the span of the taps is beyond the FFT data length).

enter image description here

$\endgroup$
  • $\begingroup$ The window method worked very well and I will try to implement the other methods you mentioned as well. Thanks for the excellent explanation. $\endgroup$ – skyguy126 Feb 28 '17 at 23:42
  • $\begingroup$ Can you gift us with an updated animation in your posting so we can see the before and after? $\endgroup$ – Dan Boschen Mar 1 '17 at 0:09
  • $\begingroup$ (If you do the other methods I think you might need more details on what I did... if interested I can paste the code in my response) $\endgroup$ – Dan Boschen Mar 1 '17 at 0:10
  • $\begingroup$ I updated the original post. As I said below this is the first time I'm working with signals in general so the results might not be as good. If you could post some example code for this segment "You can recover most of the signal by root sum-squaring the adjacent bins (or all the signal by rss'ing all the bins and displaying this at the max signal location)." it would be very helpful. $\endgroup$ – skyguy126 Mar 1 '17 at 1:30
  • 1
    $\begingroup$ I didn't do anything on that path but if you look at the graph I posted above labeled "bank of filters" you can see how the power is conserved across all the bins. If you choose any frequency and took the values in each of the 4 bins and root sum squared them $out = \sqrt{x_0^2+x_1^2+x_2^2+x_3^2}$ you would see that the result would always be 4. In a larger FFT such as yours, most of the energy would be in the adjacent bins, so you could simply root sum square over those to recover most of the signal. $\endgroup$ – Dan Boschen Mar 1 '17 at 1:39
2
$\begingroup$

You are not interpolating the magnitude peaks when they are between FFT result bin centers (basis vector frequencies). Any frequency that is not exactly integer periodic in the FFT length will be between bins, and thus split it energy among bins, which lowers the peak value in the nearest bin. Thus you see a so-called "scalloping loss".

So, you could do the opposite of damping down peaks, and instead interpolate up the peaks "hidden" between FFT result bin centers.

Or you could try applying a "flat top" window before the FFT if the FFT is long enough.

$\endgroup$
  • $\begingroup$ How would I go about interpolating the data since it is a pretty irregular curve (I assume you mean interpolate it after filling the buffer with PCM data but before the FFT). This is the first time I'm working with signal processing and audio in general. $\endgroup$ – skyguy126 Feb 26 '17 at 18:09
  • $\begingroup$ Sinc-kernel interpolation of the FFT results near the peak. $\endgroup$ – hotpaw2 Feb 26 '17 at 18:10

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.