Prev: Any good practical book on Adaptive filtering?
Next: Analyzing an "undersampled" sequence
From: Afinko on 14 Jul 2010 05:31
>If you're only concerned with amplitude/phase of a specific frequency,
>you don't need to go to the trouble of calculating a full FFT. Google
>"Goertzel algorithm." One thing to remember is that phase is relative
>to the time at which you started sampling the sinusoid; unless your
>blocks of input data to the FFT are synchronized in some way, the
>sinusoid's initial phase is going to be arbitrary.

Hi Jason,
Thanks for your reply.
It is a great idea, because I can always use so many inputs to the
"Goertzel algorithm" that I will always take the integer number of signal
periods. I used "goertzel" function in MATLAB. It works great on simulated
data. But when I used the measured data, the phase information is still
very unstable :(

The only way how I obtained gut results on measured data, that I found till
now, is this approach:
I filtered the input data with sharp bandpass IIR filter with center
frequency that I measured with laser. Then I used synchronized averaging
with reference from laser. I obtained one period of examined signal (it
looks very similar to the one sinus period). Then I simply find the maximum
of this signal and obtain a phase.
The only problem is, that this kind of IIR filter has very nonlinear phase
and difference between delay of 50 Hz and 50.2Hz is more than one period of
50 Hz signal :(

I am going to try filter the measured signal, as well as reference signal
from laser with the same IIR. In this way, the shift of examined frequency
should by the same for both signals. Hopefully it will work.

Any other ideas?
Thanks in advance.
Afinko
From: Afinko on 14 Jul 2010 05:39
>if the exact frequency, f0, is known in advance, i would just
>correlate the sampled signal to cos(2*pi*f0*n) and sin(2*pi*f0*n),
>after LPFing both separately, apply the 4 quadrant arg{} for phase
>(amplitude is from square rooting the sums of squares).

hi robert bristow-johnson,
Thanks for your reply.
The problem is, that in measured signal, there is not only this one
examined frequency component, but there is a lot of another frequency
components, with much higher amplitudes. The best solution would be to
filter this signal first, but such a sharp bandpass filter is a problem to
achieve with FIR, therefore I used IIR, as I posted earlier. The next
problem is, that the measured frequency from laser is changing for all
machine. Therefore I need to compute the filter coefficients realtime in
DSP before filter. To compute the one biquad of IIR such as with
"iirpeak.m" function in MATLAB looks quit easy.

Afinko
From: Greg Heath on 14 Jul 2010 05:40
On Jul 14, 5:16 am, "Afinko" <afinko(a)n_o_s_p_a_m.gmail.com> wrote:
> >If you reference the phase to the center of your FFT aperture
> >(by means of an pre-fftshift, or post flipping of alternate
> >bins), the evenness to oddness ratio of sinusoid can be
> >interpolated, even if it is not an integer number of signal
> >periods in the FFT aperture.  And the evenness to oddness
> >ratio is cleanly related to the atan2() phase.
>
> >The idiocy where people say you can't interpolate phase
> >is because they are using the start of the FFT aperture;
> >and a non-integer period sinusoid in discontinuous at the
> >boundaries of an FFT aperture, so trying to reference the
> >phase to a discontinuity is certainly non-pictorial and
> >non-intuitive.
>
> >Also, most windows go to zero at the boundaries of the
> >FFT aperature.  So you end up looking for the phase of
> >nothing.  The signal, after windowing, is in the middle.
>
> >So just move your phase reference to the center.
>
> Hi Ron N.,
> Thanks for your reply.
> I read your post several times, but I do not understand what exactly you
> mean by "move your phase reference to the center".
> Can you pleas post some MATLAB example how to do this?
> Or can you modify my MATLAB code and show how does it work?

For a phase reference at t = t0:

f = (Fs/N)*[0:N-1]';
X = exp(i*2*pi*f*t0).*fft(x);

Hope this helps.

Greg
From: Greg Heath on 14 Jul 2010 05:45
On Jul 13, 5:35 am, "Afinko" <afinko(a)n_o_s_p_a_m.gmail.com> wrote:
> Hi,
>
> I am trying to compute the amplitude and phase of signal with specific
> frequency. The frequency of the signal is known (measured with laser - the
> rotation machine). I tried to simulate the computation in matlab first:
>
> close all
> clear all
> % clc
>
> f = 50;                 % Frequency of signal
> fs = 64000;             % Sample frequency
> % fs = 2560;
> phase_deg = 10;         % Phase of signal [degre]
> t = 0.5;                % Time of signal
>
> x = (0:1/(fs):t-(1/fs));            % Sampled time
> phase_rad = 2*pi*phase_deg/360;     % Phase in radians
> y = cos(2*pi*f*x + phase_rad);      % Generate sampled signal
> figure;plot(x,y);                   % Plot signal
>
> % y = y + (0.1*randn(size(y)));       % Add some noise
>
> % Apply windowing function on the signal
> w = 2.*hanning(length(y));
> y = y.*w';
> figure;plot(x,y);
>
> % Compute FFT and plot abs and angle
> Y = fft(y);
> figure;plot(abs(Y));
> figure;plot(angle(Y));
>
> % Find the signal in frequency domain
> [max_num,max_pos] = max(abs(Y));
>
> % Value of signal (Amplitude)
> value_x = abs(Y(max_pos))/(length(x)/2)
>
> % Angle of signal [degre]
> angle_x = 360*(angle(Y(max_pos)))/(2*pi)
>
> I used the FFT, found the main peak (that represent the signal of
> interest), and compute the amplitude and phase of this peak.
>
> Everything works great, even a small noise don't change the value of
> amplitude and phase reasonable.
>
> However, the problem occurs when the FFT is computed from Non-Integer
> number of signal periods. What is in praxis all the time. Try to change the
> frequency of the signal to e.g. 50.4, 51.7, 53.9 etc. The classical
> approach should be to use the windowing function. I used Hanning (i tried
> lot of them) in example. Windowing will help with amplitude. But the phase
> of the signal is computed totally incorrectly, even when the windowing is
> used. I found that the windowing function do not play any (or not
> reasonable) role in phase computation.
>
> Is there any way how to compute the phase in praxis, when the FFT is
> computed from Non-Integer number of signal periods?

Replace "praxis" with "practice"?

I'm not absolutely sure of this; but maybe the
frequency interpolation obtained from time domain
zeropadding will work.

Hope this helps.

Greg
From: steve on 14 Jul 2010 11:38
On Jul 13, 5:35 am, "Afinko" <afinko(a)n_o_s_p_a_m.gmail.com> wrote:
> Hi,
>
> I am trying to compute the amplitude and phase of signal with specific
> frequency. The frequency of the signal is known (measured with laser - the
> rotation machine). I tried to simulate the computation in matlab first:
>
> close all
> clear all
> % clc
>
> f = 50;                 % Frequency of signal
> fs = 64000;             % Sample frequency
> % fs = 2560;
> phase_deg = 10;         % Phase of signal [degre]
> t = 0.5;                % Time of signal
>
> x = (0:1/(fs):t-(1/fs));            % Sampled time
> phase_rad = 2*pi*phase_deg/360;     % Phase in radians
> y = cos(2*pi*f*x + phase_rad);      % Generate sampled signal
> figure;plot(x,y);                   % Plot signal
>
> % y = y + (0.1*randn(size(y)));       % Add some noise
>
> % Apply windowing function on the signal
> w = 2.*hanning(length(y));
> y = y.*w';
> figure;plot(x,y);
>
> % Compute FFT and plot abs and angle
> Y = fft(y);
> figure;plot(abs(Y));
> figure;plot(angle(Y));
>
> % Find the signal in frequency domain
> [max_num,max_pos] = max(abs(Y));
>
> % Value of signal (Amplitude)
> value_x = abs(Y(max_pos))/(length(x)/2)
>
> % Angle of signal [degre]
> angle_x = 360*(angle(Y(max_pos)))/(2*pi)
>
> I used the FFT, found the main peak (that represent the signal of
> interest), and compute the amplitude and phase of this peak.
>
> Everything works great, even a small noise don't change the value of
> amplitude and phase reasonable.
>
> However, the problem occurs when the FFT is computed from Non-Integer
> number of signal periods. What is in praxis all the time. Try to change the
> frequency of the signal to e.g. 50.4, 51.7, 53.9 etc. The classical
> approach should be to use the windowing function. I used Hanning (i tried
> lot of them) in example. Windowing will help with amplitude. But the phase
> of the signal is computed totally incorrectly, even when the windowing is
> used. I found that the windowing function do not play any (or not
> reasonable) role in phase computation.
>
> Is there any way how to compute the phase in praxis, when the FFT is
> computed from Non-Integer number of signal periods?
>
> Thanks in advance.
> Afinko

use a single point DFT that tracks with frequency, a tracking filter
First  |  Prev  |  Next  |  Last
Pages: 1 2 3 4
Prev: Any good practical book on Adaptive filtering?
Next: Analyzing an "undersampled" sequence