How to compute phase of FFT signal in real world?
in [DSP]
|
From: Afinko on 13 Jul 2010 05:35 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
From: Jason on 13 Jul 2010 08:21 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 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. Jason
From: robert bristow-johnson on 13 Jul 2010 11:37 On Jul 13, 8:21 am, Jason <cincy...(a)gmail.com> wrote: > On Jul 13, 5:35 am, "Afinko" <afinko(a)n_o_s_p_a_m.gmail.com> wrote: > .... > > > Is there any way how to compute the phase in praxis, when the FFT is > > computed from Non-Integer number of signal periods? > > > 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." 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). r b-j
From: Ron N. on 14 Jul 2010 00:41 On Jul 13, 2: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? 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. IMHO. YMMV. -- rhn A.T nicholson d.0.t C-o-M http://www.nicholson.com/rhn/dsp.html
From: Afinko on 14 Jul 2010 05:16 >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? Thanks in advance. Afinko
|
Next
|
Last
Pages: 1 2 3 4 Prev: Any good practical book on Adaptive filtering? Next: Analyzing an "undersampled" sequence |