Prev: how to draw an arc between two points given centre and radius
Next: polyfit question
From: Heinrich Acker on 30 Mar 2010 15:23
Hello fellow Matlab users,

I am trying to find the exact frequency of an input signal (by spectral analysis). The measured portion of the signal is quite short in the sense that only samples which represent a few signal cycles are available. So I get a coarse frequency step from the fft() result. In order to get a higher resolution, I multiplied my signal s with factors of the form exp(-2*pi*i*f*t), which represent frequencies between those in the fft() output.

It looks like this:

ta = (0:numel(s)-1)/numel(s);
m = 0;
for k = 0:1/res:fmax % res=1 is fft, tried res=10 and higher
m = m+1;
A(m) = sum( s.* exp(-j*2*pi*k*ta) );
end

From playing around with this, I made two observations:

1) The method works in the sense that the maximum of the resulting spectrum shows the unknown signal frequency much better. I am not bound to the frequencies the fft() function gives. So the resolution *seems* to be much higher, just as the frequency raster that is chosen.
2) If a signal with two components is analysed like that, the resolution *seems* to be the same as with fft(), because the method can not show distinct peaks if the signal frequencies are closer than two adjacent elements of the fft() output.

What I would like to know: Is there an algorithm (interpolation-like) that can give the desired frequency from samples of the original fft() spectrum, without the additional, expensive dft-like computations in the code snipped above? I should add that I expect a dominating single sine in my application (observation 2 is not relevant for the application), well above white noise level. Is there perhaps a better method available to get a good frequency estimation from sample sets too short for dft?

Thanks a lot.

Heinrich
From: David Young on 30 Mar 2010 15:40
It's non-trivial, I'm afraid. The resolution of the DFT is what it is because intermediate frequencies convey no extra information about the signal, assuming that you're fitting sine waves to data using least squares (which is, in effect, what the DFT does). In order to get super-resolution, you need some other criterion for a good spectrum - and I know that at one time the entropy of the spectrum was the leading contender for such a criterion.

Things may have moved on, but it's probably worth looking up maximum-entropy spectral analysis and super-resolution - there may be code available too.
From: Greg Heath on 30 Mar 2010 15:40
On Mar 30, 3:23 pm, "Heinrich Acker" <firstname.lastn...(a)web.de>
wrote:
> Hello fellow Matlab users,
>
> I am trying to find the exact frequency of an input signal (by spectral analysis). The measured portion of the signal is quite short in the sense that only samples which represent a few signal cycles are available. So I get a coarse frequency step from the fft() result. In order to get a higher resolution, I multiplied my signal s with factors of the form exp(-2*pi*i*f*t), which represent frequencies between those in the fft() output.
>
> It looks like this:
>
> ta = (0:numel(s)-1)/numel(s);
> m = 0;
> for k = 0:1/res:fmax % res=1 is fft, tried res=10 and higher
> m = m+1;
> A(m) = sum( s.* exp(-j*2*pi*k*ta) );
> end
>
> From playing around with this, I made two observations:
>
> 1) The method works in the sense that the maximum of the resulting spectrum shows the unknown signal frequency much better. I am not bound to the frequencies the fft() function gives. So the resolution *seems* to be much higher, just as the frequency raster that is chosen.
> 2) If a signal with two components is analysed like that, the resolution *seems* to be the same as with fft(), because the method can not show distinct peaks if the signal frequencies are closer than two adjacent elements of the fft() output.
>
> What I would like to know: Is there an algorithm (interpolation-like) that can give the desired frequency from samples of the original fft() spectrum, without the additional, expensive dft-like computations in the code snipped above? I should add that I expect a dominating single sine in my application (observation 2 is not relevant for the application), well above white noise level. Is there perhaps a better method available to get a good frequency estimation from sample sets too short for dft?
>

1. Type in the MATLAB command window

doc goertzel
help goertzel

2.http://en.wikipedia.org/wiki/Goertzel_algorithm

3. Search this group for

goertzel

4. Extend the search to comp.dsp


Hope this helps.

Greg
From: David Young on 30 Mar 2010 15:59
Greg Heath <heath(a)alumni.brown.edu> wrote in message <0f160983-f1a3-4503-b41a-6ca4fb42d03f(a)o24g2000vbo.googlegroups.com>...
....
> 2.http://en.wikipedia.org/wiki/Goertzel_algorithm
....

Greg - are you sure Goertzel helps? The question is how to find the frequency of the sine wave when there are few cycles in the data. The Goertzel algorithm requires the frequency (or frequencies) to be specified in advance, as I understand it, and doesn't tell you any more than the DFT can.

I think what is needed is this

http://www.eos.ubc.ca/research/cdsst/Tad_home/max_ent.pdf

or, probably, more recent work that builds on this foundation. The old paper is worth looking at though, because it sets out the issue so clearly.
From: TideMan on 30 Mar 2010 17:11
On Mar 31, 8:23 am, "Heinrich Acker" <firstname.lastn...(a)web.de>
wrote:
> Hello fellow Matlab users,
>
> I am trying to find the exact frequency of an input signal (by spectral analysis). The measured portion of the signal is quite short in the sense that only samples which represent a few signal cycles are available. So I get a coarse frequency step from the fft() result. In order to get a higher resolution, I multiplied my signal s with factors of the form exp(-2*pi*i*f*t), which represent frequencies between those in the fft() output.
>
> It looks like this:
>
> ta = (0:numel(s)-1)/numel(s);
> m = 0;
> for k = 0:1/res:fmax % res=1 is fft, tried res=10 and higher
>    m = m+1;
>    A(m) = sum( s.* exp(-j*2*pi*k*ta) );
> end
>
> From playing around with this, I made two observations:
>
> 1) The method works in the sense that the maximum of the resulting spectrum shows the unknown signal frequency much better. I am not bound to the frequencies the fft() function gives. So the resolution *seems* to be much higher, just as the frequency raster that is chosen.
> 2) If a signal with two components is analysed like that, the resolution *seems* to be the same as with fft(), because the method can not show distinct peaks if the signal frequencies are closer than two adjacent elements of the fft() output.
>
> What I would like to know: Is there an algorithm (interpolation-like) that can give the desired frequency from samples of the original fft() spectrum, without the additional, expensive dft-like computations in the code snipped above? I should add that I expect a dominating single sine in my application (observation 2 is not relevant for the application), well above white noise level. Is there perhaps a better method available to get a good frequency estimation from sample sets too short for dft?
>
> Thanks a lot.
>
> Heinrich

IMHO, Fourier analysis is completely wrong in this situation.
Fourier analysis has two assumptions:
1. that the signal is stationary, which is very unlikely with a short
signal, as you describe; and
2. that the signal repeats itself ad infinitum.

From a practical point of view, unless the length of your sample is an
exact multiple of the period, there will be leakage into adjacent
frequencies and you will not get an accurate estimate of the period
(or frequency).

I'd suggest an entirely different approach called phase coherence,
described here:
Lindstrom, J; Kokko, H. and Ranta, E. 1997: Detecting periodicity in
short and noisy time series data. OIKOS, 78:2, 406-410.

It's based on chaos theory and Poincare mapping. The theory is
complicated, but the application is pretty simple. It involves just a
few lines of code to calculate the phase coherence for various subset
lengths, then you look for negative spikes in coherence which indicate
periodicity. I've found it works very well for a chaotic signal of
seiche where the period was fixed, but the amplitude and phase were
continually changing as a result of nonlinear effects (a definition of
chaos).
 |  Next  |  Last
Pages: 1 2
Prev: how to draw an arc between two points given centre and radius
Next: polyfit question