FFT Phase Plot
Brandon Helfield
I have a question about FFT phase plots that I couldn't seem to find a solution to from the other posts.
If we have a sinusoidal signal, say y=cos(2*pi*25.*t) , where t is some time vector.
Shouldn't the phase plot of the fourier transform of this signal be primarily 0 radians at 25 Hz? I understand there might be some significant phase contributions from other frequencies due to windowing (say it's top-hat windowed), but the phase plots looks crazy. Unwrapping it doesn't make sense either, because it just adds pi to all the jumps, accumulating ridiculously large phases. Having a phase of 1000 rad? That doesn't make sense.
The way I'm picturing it in my head is the phase of each of the contributions used to build up the time domain signal. Surely, a pure (windowed) cosine at f=f0 would have a phase component of 0 at f0, and most likely 0 everywhere else too.
I've also tried increasing the freq resolution. Essentially just by zero padding my time domain signal. Not only is it still not correct, but it's different from it's un-zero padded counterpart.
So confused, can anyone offer some insight?
Matt J
> Hi everyone, and thanks in advance for reading this.
>
> I have a question about FFT phase plots that I couldn't seem to find a solution to from the other posts.
>
> If we have a sinusoidal signal, say y=cos(2*pi*25.*t) , where t is some time vector.
> Shouldn't the phase plot of the fourier transform of this signal be primarily 0 radians at 25 Hz? I understand there might be some significant phase contributions from other frequencies due to windowing (say it's top-hat windowed), but the phase plots looks crazy. Unwrapping it doesn't make sense either, because it just adds pi to all the jumps, accumulating ridiculously large phases. Having a phase of 1000 rad? That doesn't make sense.
Try applying this function to the phase
atan(tan(phase))
This will bound it to the interval [-pi , pi]
Matt J
>
> I have a question about FFT phase plots that I couldn't seem to find a solution to from the other posts.
>
> If we have a sinusoidal signal, say y=cos(2*pi*25.*t) , where t is some time vector.
> Shouldn't the phase plot of the fourier transform of this signal be primarily 0 radians at 25 Hz? I understand there might be some significant phase contributions from other frequencies due to windowing (say it's top-hat windowed), but the phase plots looks crazy. Unwrapping it doesn't make sense either, because it just adds pi to all the jumps, accumulating ridiculously large phases. Having a phase of 1000 rad? That doesn't make sense.
>
> The way I'm picturing it in my head is the phase of each of the contributions used to build up the time domain signal. Surely, a pure (windowed) cosine at f=f0 would have a phase component of 0 at f0, and most likely 0 everywhere else too.
Another thing to keep in mind is that the computation of phase ignores the magnitude of the frequency component. Therefore, if there are other frequency components, no matter how small, due to numerical noise, they can have a significant contribution to the phase.
For example, this
y=cos(2*pi*25.*t) +(1e-11)* cos( 2*pi*30.*(t-t0) )
will have a phase component of 2*pi*30*t0 at 30 Hz even though the contribution of the 30 Hz frequency component is ridiculously small
Bottom line, if you have a priori knowledge of where the spectrum is supposed to be zero, you need to use it.
Greg Heath
<brandon_helfi...@hotmail.com> wrote:
> Hi everyone, and thanks in advance for reading this.
>
> I have a question about FFT phase plots that I couldn't
seem to find a solution to from the other posts.
>
> If we have a sinusoidal signal, say y=cos(2*pi*25.*t) ,
where t is some time vector.
> Shouldn't the phase plot of the fourier transform of
>this signal be primarily 0 radians at 25 Hz?
Yes
> I understand
there might be some significant phase contributions from
other frequencies due to windowing (say it's top-hat
windowed), but the phase plots looks crazy. Unwrapping
it doesn't make sense either, because it just adds pi to
all the jumps, accumulating ridiculously large phases.
Having a phase of 1000 rad? That
> doesn't make sense.
Yes, it doesn't.
That is why a "significant" phase contribution must
be defined as only occuring at significant magnitudes.
The phase is determined by the ratio of imaginary to
real parts of the transform. Therefore they will always
yield a phase in (-pi,pi] even when both are just random
noise. For example, cut and paste
eps
eps2 = eps^2
(180/pi)*atan2(eps2,-eps)
(180/pi)*atan2(eps,eps2)
(180/pi)*atan2(0,eps)
(180/pi)*atan2(-eps,eps2)
(180/pi)*atan2(-eps2,-eps)
> The way I'm picturing it in my head is the phase of
each of the contributions used to build up the time
domain signal. Surely, a pure (windowed) cosine at
f=f0 would have a phase component of 0 at f0, and most
>likely 0 everywhere else too.
No.
As I said above, even noise from floating point
error will cause the fft to have a phase in (-pi,pi].
> I've also tried increasing the freq resolution.
> Essentially just by zero padding my time domain signal.
No.
That does not increase resolution. Resolution is
predominantly determined by the length of the time
window:
df = 1/(max(t)+dt)
zero padding just INTERPOLATES at smaller frequency
intervals.
> and the Not only is it still not correct, but it's
> different from it's un-zero padded counterpart.
>
> So confused, can anyone offer some insight?
Z = fft(z)/N;
A = abs(Z);
P = angle(Z);
Am = max(max(A));
thresh = Am/100 % or whatever
mask = find(A<thresh);
Z(mask) = 0;
Hope this helps,
Greg
Greg Heath
> "Brandon Helfield" <brandon_helfi...@hotmail.com> wrote in message <hggfv5$ij...@fred.mathworks.com>...
fft phase is already bound to (-pi,pi]
Hope this helps.
Greg
Matt J
Sorry, I mean it will bound it to [-pi/2, pi/2]
Rune Allnor
wrote:
> Hi everyone, and thanks in advance for reading this.
>
> I have a question about FFT phase plots that I couldn't seem to find a solution to from the other posts.
>
> If we have a sinusoidal signal, say y=cos(2*pi*25.*t) , where t is some time vector.
> Shouldn't the phase plot of the fourier transform of this signal be primarily 0 radians at 25 Hz?
No. It should be phi = 0 + n*2*pi, where n is an
arbitrary integer.
> I understand there might be some significant phase contributions from other frequencies due to windowing (say it's top-hat windowed),
Did you use a weighting window? If so, all bets regarding
phase are off.
> but the phase plots looks crazy. Unwrapping it doesn't make sense either, because it just adds pi to all the jumps, accumulating ridiculously large phases. Having a phase of 1000 rad? That doesn't make sense.
Why not?
> The way I'm picturing it in my head is the phase of each of the contributions used to build up the time domain signal. Surely, a pure (windowed) cosine at f=f0 would have a phase component of 0 at f0, and most likely 0 everywhere else too.
>
> I've also tried increasing the freq resolution. Essentially just by zero padding my time domain signal. Not only is it still not correct, but it's different from it's un-zero padded counterpart.
>
> So confused, can anyone offer some insight?
Just do the basic trigonometry, where you keep the
2pi ambiguities in the formulas. There is nothing
more to it than that.
Rune
Matt J
> Shouldn't the phase plot of the fourier transform of this signal be primarily 0 radians at 25 Hz?
Here's perhaps a little more insight.
Your windowed cosine function is real and even, which means its spectrum is real and even as well. Moreover, the spectrum, because of the windowing, is likely to oscillate from positive to negative. For example, a rectangular window would give a spectrum consisting of two sincs centered at +/- 25 Hz.
The problem in cases like these is that for real spectra, the phase is always 0 where the spectrum is non-negative and pi where it is negative. So, you will see a phase plot that jumps discontinuously between 0 and pi wherever the spectrum oscillates from positive to negative.
Additionally, many people forget that even when the fft output is expected to be real, there will still be a small residual imaginary part due to numerical roundoff error. You can get rid of this by doing real(fft(y)). However, if you forget to do this then the phase can jump wildly and discontinuously between pi and -pi in the negative regions of the spectrum due to sign changes in the imaginary residual part.
Notice for example, how a small imaginary part can cause a great difference in the value of the phase below:
>> angle(-1 - i*1e-25)
ans =
-3.1416
>> angle(-1 + i*1e-25)
ans =
3.1416
Bottom line. If you you strip away the imaginary part, a lot of this difficulty goes away.
Matt J
Here is some code that may help illustrate what I'm saying.
It basically plots the phase spectrum of your cos(2*pi*25*t) after appropriately stripping away the imaginary part. As you will see, the phase is quite well behaved.
In addition, it shows that after doing so the phase is exactly 0 everywhere the spectrum is positive and exactly pi everywhere the spectrum is negative (see figures 3 & 4).
func=@(t) cos(2*pi*25.*t);
T=.5; %duration
N=1e3 +1; %odd number of samples
dt=T/(N); %t axis sampling interval
df=1/N/dt; % f-axis sampling interval
%set up axes
Axis= (0:N-1) -ceil((N-1)/2);
tAxis=dt*Axis; %x-axis
fAxis=df*Axis; %frequency axis
%the actual data and transform operation
y=func(tAxis);
f=(real(fftshift(fft(ifftshift(y))))); %STRIP AWAY IMAGINARY PART!!!!
%%the plot
figure(1); plot(tAxis,y)
xlabel 'Time (sec)', ylabel 'Signal Value'
figure(2); plot(fAxis, unwrap(angle(f),pi-1000*eps))
xlim([-60,60])
xlabel 'Frequency (Hz)'
ylabel 'Phase'
figure(3); plot(fAxis(f>=0), angle(f(f>=0)) )
%figure(2); plot(fAxis, (atan(tan(f))) )
xlim([-60,60])
title 'Masked to Negative Part of Spectrum'
xlabel 'Frequency (Hz)'
ylabel 'Phase'
figure(4); plot(fAxis(f<0), angle(f(f<0)))
%figure(2); plot(fAxis, (atan(tan(f))) )
xlim([-60,60])
title 'Masked to Positive Part of Spectrum'
xlabel 'Frequency (Hz)'
ylabel 'Phase (masked to negative spectrum)'
Greg Heath
> On 18 Des, 19:00, "Brandon Helfield" <brandon_helfi...@hotmail.com>
> wrote:
>
> > Hi everyone, and thanks in advance for reading this.
>
> > I have a question about FFT phase plots that I couldn't seem to find a solution to from the other posts.
>
> > If we have a sinusoidal signal, say y=cos(2*pi*25.*t) , where t is some time vector.
> > Shouldn't the phase plot of the fourier transform of this signal be primarily 0 radians at 25 Hz?
>
> No. It should be phi = 0 + n*2*pi, where n is an
> arbitrary integer.
??
angle(fft(y)) only yields values in (-pi,pi].
Hope this helps.
Greg
Brandon Helfield
Hi everyone
Thanks so much for your replies!
I'm going to try all your suggestions, and I've tried a couple.
Rune --> sure, I agree that it should be phi = 0 +2npi, but surely for a more arbitrary function (i.e not just a pure cos wave) , the phase plot should be able to give you a more appropriate looking readout, and not force you to go through to math to figure everything out.
Thanks again for everyone's help, I'm going to try your suggestions out!
Brandon
Greg Heath
wrote:
> "Matt J " <mattjacREM...@THISieee.spam> wrote in message
><hgjbe2$7l...@fred.mathworks.com>...
>
> Here is some code that may help illustrate what I'm saying.
>
> It basically plots the phase spectrum of your cos(2*pi*25*t)
> after appropriately stripping away the imaginary part. As you
> will see, the phase is quite well behaved.
>
> In addition, it shows that after doing so the phase is exactly
> 0 everywhere the spectrum is positive and exactly pi everywhere
> the spectrum is negative (see figures 3 & 4).
>
> func=@(t) cos(2*pi*25.*t);
>
> T=.5; %duration
> N=1e3 +1; %odd number of samples
> dt=T/(N); %t axis sampling interval
> df=1/N/dt; % f-axis sampling interval
>
> %set up axes
> Axis= (0:N-1) -ceil((N-1)/2);
> tAxis=dt*Axis; %x-axis
> fAxis=df*Axis; %frequency axis
1. Fs = 1/dt = 2002 is overkill for f1 = 25
2. df = 2 omits f1 = 25 from fAxis
3. max(find(fAxis<25)) = 513
4. find(fAxis==25) = Empty matrix: 0-by-1
5. min(find(fAxis>25)) = 514
> %the actual data and transform operation
> y=func(tAxis);
> %STRIP AWAY IMAGINARY PART!!!!
> f=(real(fftshift(fft(ifftshift(y)))));
'Why not use "Y" ?"f" is easily confused for frequency
> %%the plot
> figure(1); plot(tAxis,y)
> xlabel 'Time (sec)', ylabel 'Signal Value'
6. T corresponds to an ODD number of half-periods
of 25Hz
7. y(end) = y(1) = 0.0392
8. Therefore y is not windowed correctly for a
continuous periodic extension.
> figure(2); plot(fAxis, unwrap(angle(f),pi-1000*eps))
??
Why not just
plot(fAxis, angle(Z),'o')
> xlim([-60,60])
> xlabel 'Frequency (Hz)'
> ylabel 'Phase'
>
> figure(3); plot(fAxis(f>=0), angle(f(f>=0)) )
> %figure(2); plot(fAxis, (atan(tan(f))) )
>
> xlim([-60,60])
> title 'Masked to Negative Part of Spectrum'
Masked to Positive...?
> xlabel 'Frequency (Hz)'
> ylabel 'Phase'
>
> figure(4); plot(fAxis(f<0), angle(f(f<0)))
> %figure(2); plot(fAxis, (atan(tan(f))) )
>
> xlim([-60,60])
> title 'Masked to Positive Part of Spectrum'
Masked to Negative...?
> xlabel 'Frequency (Hz)'
> ylabel 'Phase (masked to negative spectrum)'
'Somewhat confusing.'
'Can this be done more simply without negative
time, fftshift and unwrapping?
What about
y = cos(2*pi*25*t) + eps*sin(2*pi*30*t);
Hope this helps.
Greg
Matt J
At the bottom, I've put a revision of the code with better labelling and which assures that f=25 Hz is sampled.
> 6. T corresponds to an ODD number of half-periods
> of 25Hz
> 7. y(end) = y(1) = 0.0392
> 8. Therefore y is not windowed correctly for a
> continuous periodic extension.
Well, I was going for something where the windowing induced more evident oscillations in the spectrum and corresponding oscillations in the phase. However, you can easily toggle W=inf in the code below to get a continuous periodic extension. In this case, the spectrum is essentially non-oscillatory, except for numerical noise, which induces very wild changes in phase. That does show a further case of numerical sensitivity of the phase, but one whos origin is not so evident from the plots.
> 'Somewhat confusing.'
> 'Can this be done more simply without negative
> time, fftshift and unwrapping?
The unwrapping was unnecessary, an inadvertent remnant of a previous draft of the code.
I don't think it's simpler, though, without negative time, fftshifting etc... Without negative time, you would have essentially a time-shifted/non-centered version of the signal that the code generates now and so would have to deal with a non-real spectrum and additional phase shift contributions. Since the OP was expecting phase=0 at 25 Hz, I assumed he was doing the same.
> What about
>
> y = cos(2*pi*25*t) + eps*sin(2*pi*30*t);
>
I tried it, and you can too by tweaking the code below. However, it doesn't produce any new effects, either for W=inf or W=25/70. For W=inf, the numerical noise alone is enough to produce wild phase fluctuations...
func=@(t) cos(2*pi*25*t);
W=20/75;
Window=@(t) (abs(t)<=W);
df=.2;% f-axis sampling interval
dt=.001; %t axis sampling interval
N=1/df/dt; %odd number of samples
%set up axes
Axis= (0:N-1) -ceil((N-1)/2);
tAxis=dt*Axis; %x-axis
fAxis=df*Axis; %frequency axis
%the actual data and transform operation
y=func(tAxis).*Window(tAxis);
Y=(real(fftshift(fft(ifftshift(y)))))*dt; %STRIP AWAY IMAGINARY PART!!!!
%%the plot
figure(1);
plot(tAxis,y)
xlabel 'Time (sec)', ylabel 'Signal Value'
xlim([-.5,.5]);
figure(2);
subplot(211)
plot(fAxis,Y)
xlabel 'Frequency (Hz)', ylabel 'Spectrum'
xlim([-60,60])
subplot(212)
plot(fAxis, angle(Y) )
xlim([-60,60])
xlabel 'Frequency (Hz)'
ylabel 'Phase'
figure(3); plot(fAxis(Y>=0), angle(Y(Y>=0)) )
%figure(2); plot(fAxis, (atan(tan(f))) )
xlim([-60,60])
title 'Masked to Positive Part of Spectrum'
xlabel 'Frequency (Hz)'
ylabel 'Phase'
figure(4); plot(fAxis(Y<0), angle(Y(Y<0)))
%figure(2); plot(fAxis, (atan(tan(f))) )
xlim([-60,60])
title 'Masked to Negative Part of Spectrum'
Brandon Helfield
Hi Matt
Thanks for the reply! But I still have a concern
If you set the imaginary part to zero, you are by definition changing the phase, setting it to 0. What about cos(... + theta) where theta is non-zero? Or, more realistically, what about a non ideal sinusoidal wave?
Matt J
>
> Thanks for the reply! But I still have a concern
>
> If you set the imaginary part to zero, you are by definition changing the phase, setting it to 0.
====================
Not everywhere. As I explained, making a spectrum real puts its phase to 0 where the spectrum is non-negative and pi elsewhere.
In the case you presented, though, there's nothing illegitimate about forcing the result to be real, because you know that the a symmetric cosine (windowed or otherwise) is theoretically supposed to be real. The fft does not produce a perfectly real result only due to machine error.
>What about cos(... + theta) where theta is non-zero?
The take-home message of what we've been talking about is that the phase function angle(x) is discontinuous on the non-negative real axis x<=0 in the complex plane.
So, if parts of your spectrum lie near the non-negative axis, you should expect numerical noise to cause discontinuous jumps in phase there. This means you have to intervene manually and strip out the numerical noise effects in these regions.
In the case of theta=0, where the spectrum is real and non-negative over large intervals, you have a particularly strong example of this case. In the case where theta~=0, the problem will show up less strongly, because the spectrum goes near the negative real axis only where abs(Y)=0 or where abs(Y) decays to 0....
>Or, more realistically, what about a non ideal sinusoidal wave?
What non-idealities, beyond what we've been discussing, are you thinking of? Windowed waves, such as we've been working with above already, are non-ideal sinusoids...
Matt J
> >What about cos(... + theta) where theta is non-zero?
>
> The take-home message of what we've been talking about is that the phase function angle(x) is discontinuous on the non-negative real axis x<=0 in the complex plane.
========
Sorry, that should be the "non-positive" real axis
strrep(EverythingIsaidAbove, 'non-negative','non-positive')
Brandon Helfield
Thanks for your help Matt.
I'm going to spend a bit more time mulling this stuff over, i'll let you know
Thanks again everyone
Brandon
Matt J
for n=0:1
%data
a=(-10:.1:-1);
RealNoise= ( 2*rand(size(a))-1 ) * 1e-25; %super-small
ImaginaryNoise= ( 2*rand(size(a))-1 ) * 1e-25 ; %super-small
ComplexBias=n*(rand+rand*i)/1000; % zero for n=0, otherwise, not small
AdditiveNoise=RealNoise+i*ImaginaryNoise;
figure(n+1);
subplot(211)
plot(a,angle(a+ComplexBias));
tt=['No Noise ']; if n, tt=[tt, ', Biased Away From Neg. Axis']; end
title(tt)
subplot(212)
plot(a,angle(a + ComplexBias + AdditiveNoise));
tt=['Noisy ']; if n, tt=[tt, ', Biased Away From Neg. Axis']; end
title(tt)
end
Greg Heath
> > Greg Heath <he...@alumni.brown.edu> wrote in message <d857e515-f3a6-43bb-813c-6342888cc...@a21g2000yqc.googlegroups.com>...
>
-----SNIP>
> > > What about
>
> > > y = cos(2*pi*25*t) + eps*sin(2*pi*30*t);
>
> > I tried it, and you can too by tweaking the code below.
However, it doesn't produce any new effects, either for W=inf
or W=25/70. For W=inf, the numerical noise alone is enough to
> > produce wild phase fluctuations...
-----SNIP
> Thanks for the reply! But I still have a concern
>
> If you set the imaginary part to zero, you are by definition
changing the phase, setting it to 0. What about cos(... + theta)
where theta is non-zero? Or, more realistically, what about a non
>ideal sinusoidal wave?
In general you should
1. plot real, imaginary, magnitude and phase to get a
more complete picture.
2. Form a mask based on an amplitude threshold.
Now Consider phase components of pi @ f1 = 25 and -pi/2
@ f2 = 30 in the following example:
z1 = cos(2*pi*f1*t+pi)
z2 = cos(2*pi*f1*t+pi) + eps*cos(2*pi*f2*t-pi/2);
(eps = 2.2204e-016)
Although differences in fft phase might only be expected
@ f = f2 = 30, the plots in Fig 2 show that differences
only occur @ f = 10 and 35. The reason is evident in Fig
1b where the magnitude of noise in the imaginary part
exceeds eps and degrades phase estimates.
The results of using amplitude threshold masking
are shown in Figs 3 and 4.
The code below yields the following premask and postmask
phase(deg) tabulations in the signal phase region.
Premask:
freq phase1 phase2 phase2-phase1
24 67.7 67.8 0.12
25 -189.0 -189.0 0.00
26 112.0 112.2 0.15
27 -152.6 -152.6 0.04
28 -168.5 -169.6 -1.13
29 -77.1 -77.2 -0.09
30 -43.0 -55.6 -12.61
31 110.0 110.1 0.18
Postmask:
freq phase1 phase2 phase2-phase1
24 0 0 0
25 -189.0 -189.0 0
26 0 0 0
27 0 0 0
28 0 0 0
29 0 0 0
30 0 0 0
31 0 0 0
Obviously, there is a 1 deg error @ f = f1 = 25
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all, close all, clc, k = 0
% Given
f1 = 25, f2 = 30
% Choose
Fs = 2.5*f2 % 75
df = 1
% Then
dt = 1/Fs % 0.0133
T = 1/df % 1
N = T/dt % 75
t = [0:dt:T-dt]';
% t = dt*(0:N-1)';
% t = linspace(0,T-dt,N)';
df = 1/T
f = 0:df:Fs-df;
% f = df*(0:N-1);
% f = linspace(0,Fs-df,N);
f([26, 31] ) % 25, 30
Fs - f([31, 26] ) % 45, 50
z1 = cos(2*pi*f1*t+pi);
z2 = z1 + eps*cos(2*pi*f2*t-pi/2);
z = [z1,z2];
Z = fft(z)/N;
for i = 1:2
X = real(Z);
Y = imag(Z);
A = abs(Z);
P = angle(Z);
k=k+1,figure(k)
subplot(2,2,1),hold on
plot(f,zeros(1,N),'k')
plot(f,X(:,1))
plot(f,X(:,2),'r')
ylabel('Real')
title('Spectra of z1(blue) and z2(red)')
subplot(2,2,2),hold on
plot(f,zeros(1,N),'k')
plot(f,eps*ones(1,N),'k--')
plot(f,-eps*ones(1,N),'k--')
plot(f,Y(:,1))
plot(f,Y(:,2),'r')
ylabel('Imaginary')
subplot(2,2,3),hold on
plot(f,A(:,1))
plot(f,A(:,2),'r')
ylabel('Amplitude')
xlabel('Frequency')
subplot(2,2,4),hold on
plot(f,zeros(1,N),'k')
plot(f,P(:,1))
plot(f,P(:,2),'r')
ylabel('Phase')
xlabel('Frequency')
[f(25:32)' (189/pi)*[P(25:32,:) ...
P(25:32,2)-P(25:32,1)]]
k=k+1,figure(k),hold on
subplot(211), hold on
plot(f,zeros(1,N),'k')
plot(f,P(:,1))
plot(f,P(:,2),'r')
ylabel('Phases')
title('Phase of z1(blue) and z2(red)')
subplot(212), hold on
plot(f,zeros(1,N),'k')
plot(f,P(:,2)-P(:,1),'g')
ylabel('Difference in Phases')
mask = find(A < max(max(A))/100);
Z(mask) = 0;
end
for j = k:-1:1
figure(j)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Hope this helps.
Greg
Greg Heath
24 64.5 64.6 0.119
25 -180.0 -180.0 0
26 106.7 106.8 0.143
27 -145.4 -145.3 0.039
28 -160.5 -161.5 -1.078
29 -73.4 -73.5 -0.083
30 -40.9 -52.9 -12.009
31 104.7 104.9 0.175
>
> Postmask:
>
> freq phase1 phase2 phase2-phase1
> 24 0 0 0
> 25 -189.0 -189.0 0
WHOOPS!
25 -180.0 -180.0 0
> 26 0 0 0
> 27 0 0 0
> 28 0 0 0
> 29 0 0 0
> 30 0 0 0
> 31 0 0 0
>
> Obviously, there is a 1 deg error @ f = f1 = 25
Ignore. There was a code typo (see below)
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> clear all, close all, clc, k = 0
-----SNIP
> [f(25:32)' (189/pi)*[P(25:32,:) ...
> P(25:32,2)-P(25:32,1)]]
WHOOPS!
[f(25:32)' (180/pi)*[P(25:32,:) ...
P(25:32,2)-P(25:32,1)]]
Hope this helps.
Greg
Matt J
> In general you should
>
> 1. plot real, imaginary, magnitude and phase to get a
> more complete picture.
> 2. Form a mask based on an amplitude threshold.
Unfortunately, amplitude masking won't be enough, in general. The angle() function is discontinuous all along the non-positive real axis. As my last example showed, If the spectrum has any purely negative value at all (large amplitude or otherwise), imaginary-valued noise can cause phase fluctuations which amplitude masking alone won't pick up.
Greg Heath
> Greg Heath <he...@alumni.brown.edu> wrote in message <585a3692-5932-4fd1-8f29-ef5f8b190...@t42g2000vba.googlegroups.com>...
> > In general you should
>
> > 1. plot real, imaginary, magnitude and phase to get a
> > more complete picture.
> > 2. Form a mask based on an amplitude threshold.
>
> Unfortunately, amplitude masking won't be enough, in general. The angle() function is discontinuous all along the non-positive real axis. As my last example showed, If the spectrum has any purely negative value at all (large amplitude or otherwise), imaginary-valued noise can cause phase fluctuations which amplitude masking alone won't pick
> up
Excellent point!
Hmmm...
How about a mask around the whole branch cut?
close all, clear all, clc, k = 0
Z0 =(-10:0.1:-1)'; % Negative spectrum
N = length(Z0);
f = (0:N-1)';
noise = 1e-25*[randn(N,1)+i*randn(N,1)]; %super-small
for n=0:1
bias = n*(rand+i*rand)/1000; % either zero or not small
Z = Z0 + bias;
A = abs(Z);
Amax = max(max(A));
thresh = Amax/100; % or whatever
mask1 = find(A < thresh);
Z(mask1) = 0;
X = real(Z);
Y = imag(Z);
mask2 = find(X<0 & abs(Y)<thresh);
Z(mask2) = X;
k=k+1,figure(k)
fftplot(f,Z)
tt=['No Noise '];
if n, tt=[tt, ', Biased Away From Neg. Axis']; end
title(tt)
Z = Z + noise;
A = abs(Z);
Amax = max(max(A));
thresh = Amax/100; % or whatever
mask1 = find(A < thresh);
Z(mask1) = 0;
X = real(Z);
Y = imag(Z);
mask2 = find(X<0 & abs(Y)<thresh);
Z(mask2) = X;
k=k+1,figure(k)
fftplot(f,Z);
tt=['Noisy '];
if n, tt=[tt, ', Biased Away From Neg. Axis']; end
title(tt)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fftplot(f,Z)
% Bare bones code for vectors
if min(size(f))~=1 | min(size(Z))~=1
error('f and Z must be a vectors')
end
N = length(Z);
X = real(Z);
Y = imag(Z);
A = abs(Z);
P = angle(Z);
subplot(2,2,1),hold on
plot(f,zeros(N,1),'k')
plot(f,X)
ylabel('Real')
title('FFT Spectrum')
subplot(2,2,2),hold on
plot(f,zeros(N,1),'k')
plot(f,Y)
ylabel('Imaginary')
subplot(2,2,3),hold on
plot(f,A)
ylabel('Amplitude')
xlabel('Frequency')
subplot(2,2,4),hold on
plot(f,zeros(N,1),'k')
plot(f,P)
ylabel('Phase')
xlabel('Frequency')
return
Hope this helps.
Greg
Matt J
> mask2 = find(X<0 & abs(Y)<thresh);
Yes, that would be the ultimate solution.
I still vaguely wonder how likely it is to have a spectrum that goes near the strictly negative real-axis, except in cases where the signal is symmetric and therefore has a purely real spectrum. In those cases, though, you would normally strip away the imaginary-valued noise via real(fft(Z)) at which point amplitude masking becomes sufficient.
Anyway, good to have a comprehensive solution, I suppose...
Matt J
> mask1 = find(A < thresh);
> mask2 = find(X<0 & abs(Y)<thresh);
On a side note, there's no need for find() here of course. Logical indexing would be sufficient....
Greg Heath
> Greg Heath <he...@alumni.brown.edu> wrote in message <2e422c51-e17f-498f-ace7-e49180467...@m38g2000yqd.googlegroups.com>...
> > mask1 = find(A < thresh);
> > mask2 = find(X<0 & abs(Y)<thresh);
>
> On a side note, there's no need for find() here of course. Logical indexing would be sufficient....
Thanks. I've been doing it longer than I care to admit.
Happy Holidays!
Greg
Brandon Helfield
Rune,
why would a weighting window mean all bets with regard to phase are off? Multiplying in the time domain, convolving in the freq. domain, with say a Gaussian envelope, shouldn't affect the phase... ?
Thanks for your help
B
Robert
Your code works well for sinusoids whose frequency falls slap bang in the middle of the FFT frequency bin e.g. bin size = 1Hz, input frequency = 5Hz, input phase = 1 radian.
But try shifting the input frequency away from the centre of the bin and the measured phase error increases dramatically, reaching a maximum half way between bins. At 5.1Hz the measured phase error is +30%, increasing about 30% for every 0.1 Hz increase in input frequency to a max of +150% at 5.5Hz.
Do you know any way round this, other than decreasing the bin size? Windowing doesn't seem to help.