MATLAB: How can i determine phase in FFT

fftifftphasezero padding

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.

Best Answer

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.

Related Solutions

MATLAB: How can i obtain different frequencies from a waveform using FFT

It doesn’t look correct to me either.

I would not use the Savitzky-Golay filter here, at least not with you original data. (I’m not certain what you want to do with it, anyway.)

If you want to get the frequency content of your signal, I would do something like this:

data = csvread('TEK00001.csv');
x=data(:,1);
y=data(:,2);
L = numel(x);                                                                           % Signal Length (Samples)
Fs = mean(diff(x));                                                                     % Sampling Frequency (Hz)
Fn = Fs/2;                                                                              % Nyquist Frequency (Hz)
nfft2=2^nextpow2(L);                                                                    % Desired Fourier Transform Length
ymc = y-mean(y);                                                                        % Subtract Mean (Eliminates d-c- Offset)
FTy = fft(ymc, nfft2)/L;                                                                % Fourier Transform (Scaled By Length Of Data Vector)
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                                                     % Frequency Vector (One-Sided Fourier Transform)
Iv = 1:numel(Fv);                                                                       % Index Vector
[pks,locidx] = findpeaks(abs(FTy(Iv))*2, 'MinPeakDistance',20, 'MinPeakHeight',2E-3);   % Peak Amplitudes & Indices
figure(2)
plot(Fv, abs(FTy(Iv))*2)
grid
xlim([0 1.5E-6])
title('FFT Magnitude')
grid minor
xlabel('Frequency [Hz]');
ylabel('Amplitude [mV]');
text(Fv(locidx), pks, sprintfc('%11.2E Hz\n%11.2E mV', [Fv(locidx)' pks]))

There are three peaks, at 1.1E-7 Hz, 3.2E-7 Hz, and 5.3E-7 Hz. (I see nothing of significance in the higher frequencies.) I used the findpeaks (link) function to isolate the most significant peak amplitudes and their frequencies, and text (link) function to label them. (The sprintfc function is undocumented. It produces a cell array of strings, very efficient for this sort of application. It otherwise works the same as sprintf.)

If you want to filter your signal, there are several options, the easiest being to use the bandpass (link) and related functions, introduced in R2018a. An appropriate high-frequency cutoff would be 1E-6 Hz.

MATLAB: 50 Hz interference and noise reduction from ECG.

hello

If you want only to remove the 50 Hz tone , use a notch filter as shown in example below. I used your data and guessed the sampling frequency from the fact that the peak was at 50 Hz.

So sampling at 10 kHz is pretty high for an ECG signal, you could easily go down by factor 10.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
signal = importdata('Sinal1.txt');
dt = 1e-4; % 0.1 milli seconds
samples = length(signal);
Fs = 1/dt; % sampling frequency (Hz)
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = Fs;    % 
Overlap = 0.75;
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
%% notch filter section %%%%%%
% H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
fc = 50;       % notch freq
wc = 2*pi*fc;
Q = 7;         % adjust Q factor for wider (low Q) / narrower (high Q) notch
                % at f = fc the filter has gain = 0
w0 = 2*pi*fc/Fs;
alpha = sin(w0)/(2*Q);
            b0 =   1;
            b1 =  -2*cos(w0);
            b2 =   1;
            a0 =   1 + alpha;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha;
            
% analog notch (for info)
num1=[1/wc^2 0 1];
den1=[1/wc^2 1/(wc*Q) 1];
% digital notch (for info)
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
freq = linspace(fc-1,fc+1,200);
[g1,p1] = bode(num1,den1,2*pi*freq);
g1db = 20*log10(g1);
[gd1,pd1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
gd1db = 20*log10(gd1);
figure(1);
plot(freq,g1db,'b',freq,gd1db,'+r');
title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');
legend('analog','digital');
xlabel('Fréquence (Hz)');
ylabel(' dB')
 
% now let's filter the signal 
signal_filtered = filtfilt(num1z,den1z,signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display : averaged FFT spectrum before / after notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq,fft_spectrum] = myfft_peak(signal, Fs, NFFT, Overlap);
sensor_spectrum_dB = 20*log10(fft_spectrum);% convert to dB scale (ref = 1)

[freq,fft_spectrum_filtered] = myfft_peak(signal_filtered, Fs, NFFT, Overlap);
sensor_spectrum_filtered_dB = 20*log10(fft_spectrum_filtered);% convert to dB scale (ref = 1)
figure(2),semilogx(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid on
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before notch filter','after notch filter');
xlabel('Frequency (Hz)');ylabel(' dB')
function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples)) = signal;
    signal = s_tmp;
    samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows
    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select);
freq_vector = (select - 1)*Fs/nfft;
end

Best Answer

Related Solutions

MATLAB: How can i obtain different frequencies from a waveform using FFT

MATLAB: 50 Hz interference and noise reduction from ECG.

Related Question