Estimating Low Frequencies (including amplitude and phase)

For the past month I've been chasing the DFT/FFT rabbit. Along the way I've been cut-n-pasting the name of other rabbits to chase later: Autocorrelation, Phase Vocoder, Parabolic Approximation, Wavelet Transform

Here's my problem. When using an FFT, the lower frequencies are always going to be missing. I actually found the same problem described here: Problems with the FFT/IFFT

Frankly it is worse than that. The energy of these frequencies are not missing at all. They are dissolved into the first several bins, making these bins semi-useless for interpolation.

For now, I simply want a general approximation of the Low Frequencies, including amplitude and phase. In FFT terms, I want bin #0.05, 0.10, 0.15, 0.20, 0.25 ...up to say around bin 3.50. Realistically, I understand that anything to small, say 0.05, would almost amount to Linear Regression. But in my mind, 0.5 or 0.75 should be doable.

I even toyed with changing k to a double (which was 'cough' "interesting"): best window for close frequency components

When Moses came down from the mountain, he did not carry tablets of stone on which it was commanded that in the equation: F = k*(df) = k*sample_rate/N, that k has to be an integer.

Anyway, I've decided it is time to try something else. So, what rabbit should I chase next? (even if it is close approximation)

PS: recommending a bigger sample is not the answer - I'm already sampling my available limit.

Test Case (I want a Blue approximation, even if lousy):

Rectangle (of course, lousy Green and so-so Red):

Gaussian 3.0 (only Red can be Interpolated, but decent):

Smaller Gaussian gives a lousy Green Interpolation and decent Red.

Edited to answer general questions (I'm sorry I did not do this to begin with):

What's the mathematical description of the data in the plot? I'm guessing Sum of Sinusoids (with noise).

What's the sample rate? All of my samples are at different rates, as in, each entire file size is used. Rather than Hertz/periods/cycles, the measurements are in deserialized data units. Most often I work with 8-byte doubles, but sometimes 4-byte integers. Therefore, with doubles, you can say that my Fs = fileBytes/8

How many samples can you capture? One each per file.

What's the application? I have no clue - it's for a client.

• Tell us more about your application, whether your signal can be described by a specific model (sum of sinusoids? harmonically related?), and what is the ultimate goal (evaluating amplitude/phase for some frequencies might not be necessary ultimately for what you want to achieve)... – pichenettes Feb 24 '13 at 18:19
• If you just want to interpolate the DFT output, you can just zero-pad the input. As pichenettes said, you assert that you need these fractional bin outputs, which may not be the case. If you can describe your ultimate goal, you may be guided toward a better approach. – Jason R Feb 24 '13 at 19:19
• Model = sum of sinusoids, which is why I also need amp/phase. And, I've tried zero-padding the time domain per stanford: link but that does not help with low freqs. – user3981 Feb 24 '13 at 20:39
• What's your signal to noise ratio? In zero noise, you can solve for a single sinewave with only 3 or 4 non-aliased points. In noise, the difference between 5% and 10% (or n+5% and n+10%) fragment of a sinewave might not be noticeable. – hotpaw2 Feb 24 '13 at 21:27
• A lot of noise. And unfortunately, closely packed peaks. I agree that unpacking a fragment will tend to be erroneous if too small. I will look into decimation for now. – user3981 Feb 25 '13 at 0:22

Updated/Revised: It sounds like you need something very powerful that goes beyond correlator/passive techniques. It sounds like you need an adaptive filter. There are various types of adaptive filtering techniques for spectral estimation, but they are often very dependent upon the application and types of input to said techniques. Could you elaborate more on how the data is collected, what you suspect is in the data, and exactly what you hope to accomplish?

If a good solution exists, I think it will require some type of adaptive filter as (I think) you are asking for something beyond the typical reach of conventional correlation/filter/FFT techniques.

A couple simple options that are generally useful but probably won't work: 1) The first one is good for interpolating a single DFT/FFT peak "in between" bins, but it relies on the fact that there is only one true peak and is concerned with finding the location of that peak.

2) If you want to calculate DFT bins at arbitrary frequencys (i.e. "bin 1.65"), then read up on the Goertzel algorithm as this will allow you to do that. Or zero-pad and then FFT. This should have the same basic effect, though, you still won't be able to resolve between bins 1.65 and 1.70 - they'll just be two points on the same gentle curve and not two distinct peaks.

• I'm alreay doing interploation, Goertzel is on the back burner (not really liking it as an option). MUSIC, very interesting, but wow, nuts trying to Google anything. "Music" Ugh? I'll keep trying though, thanks. – user3981 Mar 4 '13 at 23:59
• Okay, so no more MUSIC? As for data, it's a sum of sinusoids (refer to my picture), but with noise. As for my hope, "I want a Blue approximation". And if that is "beyond the typical reach of convention..." then help me think outside the FFT box. For instance, if I Linear regress each half of my signal, I have two lines that follow the Blue sinusoid. Stupid example, but hey I'm trying to think of "what if". Do you have anything pop in your head like this? Then let me know about it. Thanks. – user3981 Mar 5 '13 at 21:36
• After looking at papers, MUSIC is good for locating where the peaks of the spectrum are. It creates a "pseudospectrum" but there are other adaptive techniques that are probably better suited to estimate spectrum of a signal, but it often depends on the data. You should label your plots - I'm just now understanding what they mean - I think. What's the mathematical description of the data in the plot? What's the sample rate? How many samples can you capture? You also never mention the application - which is okay, but it's less info with which to understand the problem. – Dave C Mar 5 '13 at 23:52
• Okay, have a look at this paper: citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.41.2699 That should do what you want. Look at the spectral estimate near the end versus true and they specifically mention spectral phase. It's a Capon spectral estimator for real-valued signals. Here is some matlab source code of a generic capon estimator to get you started: read.pudn.com/downloads155/sourcecode/windows/comm/689621/CH5/… You may have to tailor that code based on the real-valued paper. If I test this on data, I'll try to post a clean answer. – Dave C Mar 6 '13 at 0:24
• To be clear, in the paper from the previous comment, I think you want the "APES" (Amplitude and Phase Estimate of a Sinusoid). Sounds appropriate. :-) – Dave C Mar 6 '13 at 0:32

In the frequency representation given by the FFT there are no frequencies missing and their is no binning of frequencies.

FFT decomposes a discrete signal in to a sum of sinusoids of certain frequency and phase. This representation is exact and has the same energy as the original signal. No information is lost, as you can do an inverse FFT and recover the original signal exactly. The binning is coming from the discretisation of the signal, not the FFT.

Also for a finite length signal, there are no frequencies lower than the one with the period equal to the length of the signal. This is from the mathematical definition of the spectrum. Put more simply, the DFT/FFT 'assumes' the signal is repeating. Zero padding of the signal therefore increases the period and therefore the lowers the first frequency.

If there are lower frequencies due to some prior knowledge, e.g. if we sampled a sine wave over half its period, then this information must be included.

One solution may be to roll your own spectrum and scale the frequency components with period longer than that of the signal. The assumption being that the low frequency component continues.

For example, let $f_S$ be the frequency you are interested in, with period = $1/f_s > 1$ longer than the period of the signal which is normalised to 1.

Now you convolve the signal with the sine and cosine waves at that frequency. No zero padding.

$$f_e = \sum_t^N f[t] \times cos[2\pi f_s t]$$ $$f_o = \sum_t^N f[t] \times sin[2\pi f_s t]$$

Then multiply by $1/ f_s$ to account for the assumption the signal outside of the sample continues like this, and multiply by 1/N to account for number of samples.

$$f_e' = f_e \times 1/ f_s \times 1/N$$ $$f_o' = f_o \times 1/ f_s \times 1/N$$

Amplitude and phase are then given in the usual way.

This can be done at whatever frequency graduations you desire.

Edit: Changes with regard to first comment.

Edit 2: Answer riddled with errors. Updated.

Edit 3: Add factoring for N

-- Edits welcome if I have this wrong.

• I think in computer code, not math formulas - sorry. So, I'm trying to see how your advice is different than Kevin McGee here. Is your 'To' the same as a DSP 'N', as in FFT size (e.g. 1024, 2048, etc.)? If yes, then is dividing by 'w' correct (i.e. 'To/w') or should it be flipped (i.e. 'w/To')? Also, is your 'fe' the real, and your 'fo' the imaginery parts? Thanks. – user3981 Mar 3 '13 at 17:49
• Was thinking in terms of radians, mistake though!. Have removed the $T_o$ etc, just assume for the length of the signal is '1'. $f_e$ is the even part of the signal (real) and $f_o$ is the odd part of the signal (imag). If N = 1024, t = 0:1/1024:1023/1024 ... I think it is similar to McGee. $\ast$ operation is MATLAB conv command. – geometrikal Mar 4 '13 at 11:44
• Maybe this will help dsp.stackexchange.com/questions/8088/… – geometrikal Mar 4 '13 at 22:56
• Thanks again, I'll definitely look into ZoomFFT. But to help me check my math per your formula, for a Signal=[1,2,3,4,5,6,7,8] and fs=0.7, I get Real: 57.499, Imag -18.683, is this correct? – user3981 Mar 5 '13 at 0:10
• Sorry. Using convolution was wrong. Answer not worth an upvote gaa. I used this in MATLAB now: f = [1,2,3,4,5,6,7,8]; t = linspace(0,7/8,8); fe = sum(f .* cos(2*pi*0.7*t)) * 1/0.7; fo = sum(f .* cos(2*pi*0.7*t)) * 1/0.7; – geometrikal Mar 5 '13 at 5:05

If you only want low frequencies, you can decimate (downsample) your input first and then run the same FFT. Lets say you have an input signal at 8kHz FS, but you are only interested in 100Hz and below. With a 512 point FFT you only get 15.625Hz per bin @ 8kHz FS.

Now, you can decimate (downsample) that 8kHz to 200Hz FS and then run the SAME FFT. You essentially throw away /LPF everything >= 100Hz. Now, each bin contains 200Hz / 512 or 0.390625Hz per bin.

You are essentially running a much longer FFT. There is extra processing, though, for the downsample, but no extra memory required for the FFT, which is probably where you run into the "max length" wall, right?

Edit: Found a link to a relevant answer: Downsampling lowpass filter for audio: FIR or IIR?

• Correct. This is what I need - the LF's to be 'per bin' a.k.a. the Hz per bin to be decimals. So, I'll start chasing the decimate rabbit this week. – user3981 Feb 25 '13 at 0:12
• Take a look at the polyphase filter link in the link that I posted. That is a pretty low-overhead way of downsampling by an integer ratio. For example, if you downsample by a factor of 8 (new FS = old FS/8) and you want to use a 256 tap FIR filter, you actually only need to use an average of 256/8 = 32 MAC operations per output sample as you only run the FIR for your final output samples, not the ones that you are throwing away during decimation. – c.fogelklou Feb 25 '13 at 7:46
• Note that decimating before your DFT won't give you any more frequency resolution; it will only give you bins with a finer spacing. Due to spectral leakage, in order to improve your ability to resolve closely-spaced sinusoids using a DFT, you need to lengthen your observation time, which you said is off the table. – Jason R Feb 25 '13 at 14:05
• @Jason R, he said that his sample length couldn't increase. I assume that he means his FFT size. If he means FFT size, then with the same size FFT, with decimation, he is effectively lengthening the observation time. He will get more resolution with the decimated signal because the FFT (of the same size) will be sampling a much lower FS (but much higher observation time.) – c.fogelklou Feb 25 '13 at 14:18
• Correct again - my DFT size is the whole signal file, my FFT is the first pow of 2 under DFT size (unless I zero-pad, then FSx5.7 to fit next pow 2). But then again, let's keep our minds open - I've mentioned that I do not have to continue chasing the DFT/FFT. So, Jason, what would you chase? As I said, I'm just looking for a general approximation here, even if it's a LF sum that I can use as a momentary mimic in time. (For now, I'll chase decimation). – user3981 Feb 25 '13 at 15:54