How can i determine phase in FFT?

Hi everyone, right now im trying to calculate signal phases using angle(x) from FFT Function im Matlab. Noted that i've coded the program like below :

%%Plotting Grafik
%create a time vector 't', containing integers from 1 to n(summary of data)

count= length(data);

Ts=mean(diff(times1));
Fs=1/Ts;

NFFT=2^nextpow2(count);
y=fft(data,NFFT)/count;

f=Fs/2*linspace(0,1,NFFT/2+1);

figure(2)
plot(f,2*abs(y(1:NFFT/2+1)));
grid on
title('Single-Sided Amplitude Spectrum of y(t)');
xlabel('Frequency(Hz)');
ylabel('|Y(f)|');
phase=angle(y);

But i'm not really sure with the phase function. Does anyone ever has the same experience?

Thanks before for your answer.

[Dan Bullard]

I use this constantly in Excel, it's atan2(real, imaginary)

https://gcc.gnu.org/onlinedocs/gfortran/ATAN2.html

With this I proved Bullard Laws of Harmonics #5,

http://www.danbullard.com/dan/articles/why_law_5_is_so_important/why_law_5_is_so_important.html

[Star Strider]

That’s backwards. See the Wikipedia article on atan2 for clarification.

[Bruno Luong]

@Dan Bullard "it's atan2(real, imaginary)"

That formula is wrong the right formula is atan2(imaginary, real) and this is exactly what angle function does.

So reverse the two atan2 arguments, then you will fix the "backwardness" of your formula.

[Dan Bullard]

arctan2 is different in Fortran and Excel, forgive me, I use Excel and didn't know about the angle() function.

[WH Wang]

it's atan2(imaginary,real) in matlab, or using phase() function of matlab.

[Star Strider]

There is no phase function that I can find in the documentation, however there is the angle function that returns the phase angle (in radians) of a complex argument.

Answer 1

[Star Strider]

Your code appears to me to be correct, including the ‘phase’ assignment. Note that the units of ‘phase’ are in radians. Convert them to degrees (if desired) by multiplying them by 180/pi.

I’m not certain what you want, so also consider the unwrap function. Use it on the radian units, and convert to degrees later if necessary.

EDIT — If you want to plot it, the easiest way is to use the subplot function. Your plotting code changes to:

phase=angle(y);

figure(2)
subplot(2,1,1)
plot(f,2*abs(y(1:NFFT/2+1)));
grid on
title('Single-Sided Spectrum of y(t)');
ylabel('|Y(f)|');
subplot(2,1,2)
plot(f,phase)
grid on
xlabel('Frequency(Hz)');
ylabel('Phase (rad)')

I did not test that, but it should work.

[Big dream]

Hi Star rider,

Thx for your help. There was a problem with plot funct, because 'phase' has array 2 times more than 'f'. Thats why i tried this code:

%Try to plot Phase FFT

phase=angle(y(1:NFFT/2)); figure(3) plot(f,phase); grid on xlabel('Frequency (Hz)') ylabel('Phase(rad)')

It still correct, right?

[Star Strider]

Yes.

I would do something similar to:

plot(f,phase(1:length(f)))

or:

plot(f,phase(1:numel(f))))

The ‘phase’ vector has to have the same length as the frequency vector.

[Big dream]

Hi.

Yes, the result is the same. Thanks by the way,this script is really useful to find freq and phase of a signal. And i've question, for the FFT function, i got a function for

y=fft(data,NFFT)/count;

instead of y=fft(data,NFFT); because i thought, if i want to convert a signal wave to freq domain, then i have to divide 'y' with length of sample data.

Or did actually the calculation 'y=fft(data,NFFT);' already the best answer? also do you know any idea how to use IFFT, to get my signal back? I've tried

yinv=ifft(data,NFFT)

but i got the wrong values.

[Star Strider]

My pleasure.

This is correct:

y = fft(data,NFFT)/count;

You have to normalise the result by the length of the original data (the ‘energy’ in the original signal) to get the correct amplitudes for the double-sided Fourier transform. (To get the correct amplitudes for the single-sided Fourier transform, you coded it correctly in multiplying it by 2.)

To get the inverse transform, take the inverse of the transformed signal:

yinv = ifft(y);

I’ve never actually done that experiment before now, so I coded it to see the result. It was as I expected, the only challenge being to define correctly the new time vector for ‘yinv’ (although I named my variables differently so as not to confuse them with your code):

Ts = 0.01;
Fs = 1/Ts;
Fn = Fs/2;
t = linspace(0, 1, 1500)*Ts;
s = sin(2*pi*t*5000) .* cos(2*pi*t*150);                        % Create Original Data

NFFT = 2^nextpow2(length(s));
FTs = fft(s, NFFT);                                             % Fourier Transform Of Zero-Padded Original Data
IFTs = ifft(FTs);                                               % Inverse Fourier Transform Of Zero-Padded Original Data
t2 = linspace(0, 1, NFFT)*Ts*NFFT/length(t);                    % Inverse Fourier Transform Sampling Interval = Ts*NFFT/length(t)

figure(1)
subplot(2,1,1)
plot(t,s)
title('Original Signal')
grid
subplot(2,1,2)
plot(t2,IFTs)
grid
title('IFFT Of FFT Of Zero-Padded Original Signal')
xlabel('Time (sec)')

Note that the times for both signals are the same, with the zero-padding clearly apparent in the second signal.

I learned something from this, so I thank you!

[Big dream]

you're welcome and Thanks again.

[Star Strider]

My pleasure.

[Big dream]

Hallo Star Rider, could i ask you personally? I've question regarding this topic. https://de.mathworks.com/matlabcentral/answers/308872-how-can-i-determine-phase-in-fft

I've read your previous comment about filtering in FFT Function https://de.mathworks.com/matlabcentral/answers/260277-how-to-do-digital-filtering-in-matlab-with-a-specified-cut-off-frequency

that's why in my programm i tried to filter my signal in time domain with a coefficient's filter from bandpass filter. But, i think its result not the correct answer. Because after applied the filter, i tried to implement the filtered signal to Frequency domain, its magnitude response go down from 20(before filter) until 0.8(after filter).

I had only using function 'result=filtfilt(coefficientsfilter,mysignal)'. The coefficient's filter comes from the function designfilt Bandpass filter.

1. Should i actually do something before i apply the filter? 2. Actually are the result from FFT function and Freqz function the same? because i've read in MATLAB page, they're almost similar.

Thanks before

[Star Strider]

I do not completely understand the problem. It seems that you are filtering the signal correctly.

Depending on the frequency content of your signal, the result of dividing the Fourier transform of the filtered signal by the Fourier transform of the unfiltered signal should closely approximate the result calculated by the freqz function.

If you are looking at your filtered signal in the time domain, the amplitude will probably decrease, depending on the frequency content of your original signal and the characteristics of your filter. The reason is that you are filtering out some of the frequencies and therefore some of the energy from the original signal.

[Big dream]

Star Rider,That's why. Because when i plot the timesamples(in Frequency and Time domain) after the filtering process, the result is worse than without filter.

Could i ask for your help, to take a look at my program below, i've also attached my timesamples. In this programm i want to remove the noise below 100kHz and above 400kHz(Bandpass 100-400kHz). I've chosen this range because my signal information has range from 100-400KHz:

Fs=4e6;

Ts=1/Fs;

count=length(timesamples);

coeffilter=designfilt('bandpassfir',... 'StopbandFrequency1', 1e5, 'StopbandFrequency2', 4e5,... 'PassbandFrequency1', 1.5e5, 'PassbandFrequency2' , 2.5e5,... 'StopbandAttenuation1', 40,'PassbandRipple',3,... 'StopbandAttenuation2', 40, 'SampleRate', 4e6);

result=filtfilt(coeffilter,timesamples);

%%end of programm

Thank you very much.

[Star Strider]

If you want your passband to be 100 kHz to 400 kHz, you have to give the correct information to the designfilt function. You set your passband at 150 kHz to 250 kHz.

Run this line to see your filter characteristics:

freqz(coeffilter, 4096, Fs)

This may be closer to the filter you want:

Fs = 4e6;
coeffilter=designfilt('bandpassfir',...
             'StopbandFrequency1', 7.5e4, 'StopbandFrequency2', 4.25e5,...
             'PassbandFrequency1', 1.0e5, 'PassbandFrequency2' , 4.0e5,...
             'StopbandAttenuation1', 40,'PassbandRipple',3,...
              'StopbandAttenuation2', 40, 'SampleRate', 4e6);

freqz(coeffilter, 4096, Fs)

[Big dream]

Thank you very much ^^. It more better, than before.

But, could you please tell me for what reason the value stopbandfrequency1 & 2 are 7.5e4 &4.25e5, do you have maybe some reference.?

[Star Strider]

My pleasure.

No reference, just my own experience. I chose those to provide a steep rolloff at the band edges. Experiment with them to see the result in the freqz plot and in your signal.

It is best that they be symmetrically placed (the same distance) from your lower and upper passband frequencies. There are no strict requirements, providing they produce a stable filter and the result you want.

[Big dream]

thanks again for your advice.^^

[Star Strider]

My pleasure.

Answer 2

[karinkan]

I am a bit confusing that the following two equations which one is correct? eq.1 y=fft(data,NFFT)/count; eq.2 y=fft(data,NFFT)/NFFT;

[Star Strider]

The first one.