Hello,
There is something I don't understand:
the discrete Fourier transform of an even real function
is expected to be even and real too, and then the phase of this dft
should be equal to 0°.
So, the phase of the dft of a gaussian function (which is also a gaussian
function), should be equal to 0°.
I wrote this on Matlab:
x = -0.6:0.001:0.6;
testx = exp(-((50*x).^2));
plot(1000*x,testx)
title('testx')
xlabel('x (millimeters)')
tftestxinit = fft(testx);
tftestx = fftshift(tftestxinit);
angletftestx = angle(tftestx);
f = 1000*(-600:600)/1201;
figure;plot(f,angletftestx)
title('angular frequency content of testx')
xlabel('f (1/m)')
ylabel('phase of the dft')
and I see a strange figure..
The phase is not equal to 0° !!?
Can somebody help me?
Thanks
mariauld
phase of fft on matlab
Started by ●March 15, 2008
Reply by ●March 15, 20082008-03-15
Hi,
The trouble comes when computing angle of complex number when
both real and imaginary part is close to zero. for example
c = a + jb
if a = 1e-17 and b = 1e-18, both of them are very close to zero
but doing atan(b/a) will not yield correct answer and will create
lot of fluctuations in the phase response. The way to clean this
up is to force all values of a and b which are lesser then certain
threshold to go to zero. This will help cleanup the phase response
of unwanted fluctuations. After this what you get is the desired
phase response.
I have modified your code. please check that it is giving the desired
result. Not sure about phase in the middle part (maybe it is expected),
but sure it has cleaned up unwanted phase jumps.
clear all;
close all;
x = -0.6:0.001:0.6;
testx = exp(-((50*x).^2));
plot(1000*x,testx)
title('testx')
xlabel('x (millimeters)')
tftestxinit = fft(testx);
tftestx = fftshift(tftestxinit);
%---------------------------- code modification start
hz = tftestx;
hzr = real(tftestx);
hzi = imag(tftestx);
index = find(abs(hzr) < 1e-10);
hzr(index) = 0;
index = find(abs(hzi) < 1e-10);
hzi(index) = 0;
tftestm = complex(hzr, hzi);
%---------------------------------- code modification ended
angletftestx = angle(tftestm);
f = 1000*(-600:600)/1201;
figure;plot(f,angletftestx)
title('angular frequency content of testx')
xlabel('f (1/m)')
ylabel('phase of the dft')
Reply by ●March 15, 20082008-03-15
Hello Bharat, Is problem the computation of the phase compute function, or that signals when they are so small may be mostley noise, so phase number is not meaningful? If you always set small numbers to zero and plug into phase compute function you always get same answer, which is not meaningful, and if caring about continuous phase, it will not be. Depends on application, if continuous phase desired, may want to interpolate between surrounding phase values that you trust. Regards, Kamar Ruptan DSP Guru bharat pathak wrote:> Hi, > > The trouble comes when computing angle of complex number when > both real and imaginary part is close to zero. for example > > c = a + jb > > if a = 1e-17 and b = 1e-18, both of them are very close to zero > but doing atan(b/a) will not yield correct answer and will create > lot of fluctuations in the phase response. The way to clean this > up is to force all values of a and b which are lesser then certain > threshold to go to zero. This will help cleanup the phase response > of unwanted fluctuations. After this what you get is the desired > phase response. > > I have modified your code. please check that it is giving the desired > result. Not sure about phase in the middle part (maybe it is expected), > but sure it has cleaned up unwanted phase jumps. > > clear all; > close all; > > x = -0.6:0.001:0.6; > testx = exp(-((50*x).^2)); > plot(1000*x,testx) > title('testx') > xlabel('x (millimeters)') > > tftestxinit = fft(testx); > tftestx = fftshift(tftestxinit); > > %---------------------------- code modification start > > hz = tftestx; > > hzr = real(tftestx); > hzi = imag(tftestx); > index = find(abs(hzr) < 1e-10); > hzr(index) = 0; > index = find(abs(hzi) < 1e-10); > hzi(index) = 0; > tftestm = complex(hzr, hzi); > > %---------------------------------- code modification ended > > angletftestx = angle(tftestm); > > f = 1000*(-600:600)/1201; > figure;plot(f,angletftestx) > title('angular frequency content of testx') > xlabel('f (1/m)') > ylabel('phase of the dft')
Reply by ●March 15, 20082008-03-15
Hello Bharat, Is problem the computation of the phase compute function, or that signals when they are so small may be mostley noise, so phase number is not meaningful? If you always set small numbers to zero and plug into phase compute function you always get same answer, which is not meaningful, and if caring about continuous phase, it will not be. Depends on application, if continuous phase desired, may want to interpolate between surrounding phase values that you trust. Regards, Kamar Ruptan DSP Guru bharat pathak wrote:> Hi, > > The trouble comes when computing angle of complex number when > both real and imaginary part is close to zero. for example > > c = a + jb > > if a = 1e-17 and b = 1e-18, both of them are very close to zero > but doing atan(b/a) will not yield correct answer and will create > lot of fluctuations in the phase response. The way to clean this > up is to force all values of a and b which are lesser then certain > threshold to go to zero. This will help cleanup the phase response > of unwanted fluctuations. After this what you get is the desired > phase response. > > I have modified your code. please check that it is giving the desired > result. Not sure about phase in the middle part (maybe it is expected), > but sure it has cleaned up unwanted phase jumps. > > clear all; > close all; > > x = -0.6:0.001:0.6; > testx = exp(-((50*x).^2)); > plot(1000*x,testx) > title('testx') > xlabel('x (millimeters)') > > tftestxinit = fft(testx); > tftestx = fftshift(tftestxinit); > > %---------------------------- code modification start > > hz = tftestx; > > hzr = real(tftestx); > hzi = imag(tftestx); > index = find(abs(hzr) < 1e-10); > hzr(index) = 0; > index = find(abs(hzi) < 1e-10); > hzi(index) = 0; > tftestm = complex(hzr, hzi); > > %---------------------------------- code modification ended > > angletftestx = angle(tftestm); > > f = 1000*(-600:600)/1201; > figure;plot(f,angletftestx) > title('angular frequency content of testx') > xlabel('f (1/m)') > ylabel('phase of the dft')
Reply by ●March 15, 20082008-03-15
Welcome Kamar,
Good to see you post with very few mistakes. Also that
you have retained the title of "DSP Guru".
What does PHD2BE mean? Are you pursuing PHD? If so what
is the topic of your research?
Regards
Bharat
>Hello Bharat,
>
>Is problem the computation of the phase compute function, or that
>signals when they are so small may be mostley noise, so phase number
>is not meaningful?
>
>If you always set small numbers to zero and plug into phase compute
>function you always get same answer, which is not meaningful, and if
>caring about continuous phase, it will not be.
>
>Depends on application, if continuous phase desired, may want to
>interpolate between surrounding phase values that you trust.
>
>Regards,
>
>Kamar Ruptan
>DSP Guru
>
>bharat pathak wrote:
>> Hi,
>>
>> The trouble comes when computing angle of complex number when
>> both real and imaginary part is close to zero. for example
>>
>> c = a + jb
>>
>> if a = 1e-17 and b = 1e-18, both of them are very close to zero
>> but doing atan(b/a) will not yield correct answer and will create
>> lot of fluctuations in the phase response. The way to clean this
>> up is to force all values of a and b which are lesser then certain
>> threshold to go to zero. This will help cleanup the phase response
>> of unwanted fluctuations. After this what you get is the desired
>> phase response.
>>
>> I have modified your code. please check that it is giving the
desired
>> result. Not sure about phase in the middle part (maybe it is
expected),
>> but sure it has cleaned up unwanted phase jumps.
>>
>> clear all;
>> close all;
>>
>> x = -0.6:0.001:0.6;
>> testx = exp(-((50*x).^2));
>> plot(1000*x,testx)
>> title('testx')
>> xlabel('x (millimeters)')
>>
>> tftestxinit = fft(testx);
>> tftestx = fftshift(tftestxinit);
>>
>> %---------------------------- code modification start
>>
>> hz = tftestx;
>>
>> hzr = real(tftestx);
>> hzi = imag(tftestx);
>> index = find(abs(hzr) < 1e-10);
>> hzr(index) = 0;
>> index = find(abs(hzi) < 1e-10);
>> hzi(index) = 0;
>> tftestm = complex(hzr, hzi);
>>
>> %---------------------------------- code modification ended
>>
>> angletftestx = angle(tftestm);
>>
>> f = 1000*(-600:600)/1201;
>> figure;plot(f,angletftestx)
>> title('angular frequency content of testx')
>> xlabel('f (1/m)')
>> ylabel('phase of the dft')
>
Reply by ●March 15, 20082008-03-15
On Mar 15, 2:06�pm, "bharat pathak" <bha...@arithos.com> wrote:> Welcome Kamar, > > � � Good to see you post with very few mistakes. Also that > � � you have retained the title of "DSP Guru". > > � � What does PHD2BE mean? Are you pursuing PHD? If so what > � � is the topic of your research? > > Regards > Bharat > > > > > > >Hello Bharat, > > >Is problem the computation of the phase compute function, or that > >signals when they are so small may be mostley noise, so phase number > >is not meaningful? > > >If you always set small numbers to zero and plug into phase compute > >function you always get same answer, which is not meaningful, and if > >caring about continuous phase, it will not be. > > >Depends on application, if continuous phase desired, may want to > >interpolate between surrounding phase values that you trust. > > >Regards, > > >Kamar Ruptan > >DSP Guru > > >bharat pathak wrote: > >> Hi, > > >> � The trouble comes when computing angle of complex number when > >> � both real and imaginary part is close to zero. for example > > >> � c = a + jb > > >> � if a = 1e-17 and b = 1e-18, both of them are very close to zero > >> � but doing atan(b/a) will not yield correct answer and will create > >> � lot of fluctuations in the phase response. The way to clean this > >> � up is to force all values of a and b which are lesser then certain > >> � threshold to go to zero. This will help cleanup the phase response > >> � of unwanted fluctuations. After this what you get is the desired > >> � phase response. > > >> � I have modified your code. please check that it is giving the > desired > >> � result. Not sure about phase in the middle part (maybe it is > expected), > >> � but sure it has cleaned up unwanted phase jumps. > > >> clear all; > >> close all; > > >> x = -0.6:0.001:0.6; > >> testx = exp(-((50*x).^2)); > >> plot(1000*x,testx) > >> title('testx') > >> xlabel('x (millimeters)') > > >> tftestxinit = fft(testx); > >> tftestx = fftshift(tftestxinit); > > >> %---------------------------- code modification start > > >> hz � � = tftestx; > > >> hzr � = real(tftestx); > >> hzi � = imag(tftestx); > >> index = find(abs(hzr) < 1e-10); > >> hzr(index) = 0; > >> index = find(abs(hzi) < 1e-10); > >> hzi(index) = 0; > >> tftestm � = complex(hzr, hzi); > > >> %---------------------------------- code modification ended > > >> angletftestx = angle(tftestm); > > >> f = 1000*(-600:600)/1201; > >> figure;plot(f,angletftestx) > >> title('angular frequency content of testx') > >> xlabel('f (1/m)') > >> ylabel('phase of the dft')- Hide quoted text - > > - Show quoted text -Hello Bharat, Topic I not want to open hear to feed the sharks. Was my response post useful to you? Regards, Kamar Ruptan DSP Guru
Reply by ●March 16, 20082008-03-16
Hello,
check this code:
%%%%%%%%%%%%%%%%%%
clear all
close all
n = [-3:0.001:3]';
x = exp(-(n.^2));
plot(n,x)
title('x')
X = fft(x);
X_mag = abs(X);
X_phase = angle(X);
figure;
plot(n,X_phase);
title('phase response')
%%%%%%%%%%%%%%%%%%%%%
this is pretty much the same thing as the initial code, except that i have
changed the variance of the signal and i have used somehow more readable
names of variables. Now, run the code. You will see that the phase is
almost perfectly linear. I have not applied none of the before-mentioned
suggestions. Now, change the number of samples, not by increasing the
sampling frequency (1000Hz) but by increasing the interval from [-3 3] to
[-30 30]. The phase response is far from linear.
Who can explain that?
Manolis
Reply by ●March 17, 20082008-03-17
On Mar 15, 7:36�am, "mariauld" <marius.pel...@institutoptique.fr> wrote:> Hello, > > There is something I don't understand: > the discrete Fourier transform of an even real function > is expected to be even and real too, and then the phase of this dft > should be equal to 0�. > So, the phase of the dft of a gaussian function (which is also a gaussian > function), should be equal to 0�. > > I wrote this on Matlab: > > x = -0.6:0.001:0.6; > testx = exp(-((50*x).^2)); > plot(1000*x,testx) > title('testx') > xlabel('x (millimeters)') > > tftestxinit = fft(testx); > tftestx = fftshift(tftestxinit); > > angletftestx = angle(tftestx); > > f = 1000*(-600:600)/1201; > figure;plot(f,angletftestx) > title('angular frequency content of testx') > xlabel('f (1/m)') > ylabel('phase of the dft') > > and I see a strange figure.. > The phase is not equal to 0� !!? > > Can somebody help me? > Thanks > > mariauldmariauld, Actually your function is not a real even function when you input it into the fft. I was defined as a real even function symmetric about x=0, but you shifted it in x by 0.6 to make the x=-0.6 sample align with n=0 for the fft. Then there is no reason to expect the phase out of the fft to be 0. If you had put the samples for 0 to 0.6 in the first part of the FFT input vector, and the samples for -0.6+0.001 to 0-0.001 in the remaining part of the input vector then the function would be real and even as far as the FFT was concerned. Note that this uses one less sample than you generated (assuming samples were generated for both -0.6 and 0.6). Dirk
Reply by ●March 17, 20082008-03-17
On Mar 17, 11:40�am, dbell <bellda2...@cox.net> wrote:> On Mar 15, 7:36�am, "mariauld" <marius.pel...@institutoptique.fr> > wrote: > > > > > > > Hello, > > > There is something I don't understand: > > the discrete Fourier transform of an even real function > > is expected to be even and real too, and then the phase of this dft > > should be equal to 0�. > > So, the phase of the dft of a gaussian function (which is also a gaussian > > function), should be equal to 0�. > > > I wrote this on Matlab: > > > x = -0.6:0.001:0.6; > > testx = exp(-((50*x).^2)); > > plot(1000*x,testx) > > title('testx') > > xlabel('x (millimeters)') > > > tftestxinit = fft(testx); > > tftestx = fftshift(tftestxinit); > > > angletftestx = angle(tftestx); > > > f = 1000*(-600:600)/1201; > > figure;plot(f,angletftestx) > > title('angular frequency content of testx') > > xlabel('f (1/m)') > > ylabel('phase of the dft') > > > and I see a strange figure.. > > The phase is not equal to 0� !!? > > > Can somebody help me? > > Thanks > > > mariauld > > mariauld, > > Actually your function is not a real even function when you input it > into the fft. �I was defined as a real even function symmetric about > x=0, but you shifted it in x by 0.6 to make the x=-0.6 sample align > with n=0 for the fft. �Then there is no reason to expect the phase out > of the fft to be 0. > > If you had put the samples for 0 to 0.6 in the first part of the FFT > input vector, and the samples for -0.6+0.001 to 0-0.001 in the > remaining part of the input vector then the function would be real and > even as far as the FFT was concerned. Note that this uses one less > sample than you generated (assuming samples were generated for both > -0.6 and 0.6). > > Dirk- Hide quoted text - > > - Show quoted text -Ignore my last post. I will post more tonight re this when I have more time. Dirk
Reply by ●March 17, 20082008-03-17
Hi,
first I'd like to thank you all, Bharat, Kamar, Manolis, Dirk, for your
useful help.
Bharat, your idea of cleaning up the numbers close to zero
helped me to avoid noise effects.
I associated this with the idea of Dirk, using fftshift on vector x.
This leads me to good results:
x1 = -0.6:0.001:0.6;
x = ifftshift(x1);
testx = exp(-((50*x).^2));
tftestxinit = fft(testx);
tftestx = fftshift(tftestxinit);
%---------------------------- code modification start
hz = tftestx;
hzr = real(tftestx);
hzi = imag(tftestx);
index = find(abs(hzr) < 1e-10);
hzr(index) = 0;
index = find(abs(hzi) < 1e-10);
hzi(index) = 0;
tftestm = complex(hzr, hzi);
%---------------------------------- code modification ended
abstftestx = abs(tftestm);
angletftestx = angle(tftestm);
f = 1000*(-600:600)/1201;
figure;plot(f,abstftestx)
title('absolute frequency content of testx')
xlabel('f (1/m)')
ylabel('magnitude of the dft')
figure;plot(f,angletftestx)
title('angular frequency content of testx')
xlabel('f (1/m)')
ylabel('phase of the dft')
Thank you all again
