Finding peak frequency in a complex signal using Matlab

1k views Asked by At

I'm trying to find the peak frequency for two signals 'CA1' and 'PFC', within a specified range (25-140Hz).

In Matlab, so far I have plotted an FFT for each of these signals (see pictures below). These FFTs suggest that the peak frequency between 25-140Hz is different for each signal, but I would like to quantify this (e.g. CA1 peaks at 80Hz, whereas PFC peaks at 55Hz). However, I think the FFT is not smooth enough, so when I try and extract the peak frequencies it doesn't make sense as my code pulls out loads of values. I was only expecting a few values - one each time the FFT peaks (around 2Hz, 5Hz and ~60Hz).

I want to know, between 25-140Hz, what is the peak frequency in 'CA1' compared with 'PFC'. 'CA1' and 'PFC' are both 152401 x 7 matrices of EEG data, recorded from 7 separate individuals. I want the MEAN peak frequency for each data set (i.e. averaged across the 7 test subjects for CA1 and PFC).

My code so far (based on Matlab help files and code I've scrabbled together online):

Fs = 508;

%notch filter
[b50,a50] = iirnotch(50/(Fs/2), (50/(Fs/2))/70);
CA1 = filtfilt(b50,a50,CA1);
PFC = filtfilt(b50,a50,PFC); 

%FFT
L = length(CA1); 

NFFT = 2^nextpow2(L);

%FFT for each of the 7 subjects
for i = 1:size(CA1,2);

    CA1_FFT(:,i) = fft(CA1(:,i),NFFT)/L;

    PFC_FFT(:,i) = fft(PFC(:,i),NFFT)/L;

end

%Average FFT across all 7 subjects - CA1
Mean_CA1_FFT = mean(CA1_FFT,2);
% Mean_CA1_FFT_abs = 2*abs(Mean_CA1_FFT(1:NFFT/2+1));

%Average FFT across all 7 subjects - PFC
Mean_PFC_FFT = mean(PFC_FFT,2);
% Mean_PFC_FFT_abs = 2*abs(Mean_PFC_FFT(1:NFFT/2+1));

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

%LEFT HAND SIDE FIGURE
plot(f,2*abs(Mean_CA1_FFT(1:NFFT/2+1)),'r');
set(gca,'ylim', [0 2]);
set(gca,'xlim', [0 200]);

[C,cInd] = sort(2*abs(Mean_CA1_FFT(1:NFFT/2+1)));
CFloor = 0.1; %CFloor is the minimum amplitude value (ignore small values)
Amplitudes_CA1 = C(C>=CFloor); %find all amplitudes above the CFloor
Frequencies_CA1 = f(cInd(1+end-numel(Amplitudes_CA1):end)); %frequency of the peaks

%RIGHT HAND SIDE FIGURE
figure;plot(f,2*abs(Mean_PFC_FFT(1:NFFT/2+1)),'r');
set(gca,'ylim', [0 2]);
set(gca,'xlim', [0 200]);

[P,pInd] = sort(2*abs(Mean_PFC_FFT(1:NFFT/2+1)));
PFloor = 0.1; %PFloor is the minimum amplitude value (ignore small values)
Amplitudes_PFC = P(P>=PFloor); %find all amplitudes above the PFloor
Frequencies_PFC = f(pInd(1+end-numel(Amplitudes_PFC):end)); %frequency of the peaks

enter image description here

Please help!! How do I calculate the 'major' peak frequencies from an FFT, and ignore all the 'minor' peaks (because the FFT is not smoothed).

1

There are 1 answers

0
Tim B On

FFTs assume that the signal has no trend (this is called a stationary signal), if it does then this will give a dominant frequency component at 0Hz as you have here. Try using the MATLAB function detrend, you may find this solves your problem.

Something along the lines of:

x = x - mean(x)
y = detrend(x, 'constant')