How can I visualize by frequency for frequency range 0 - 70.87Hz? (EEG_1=3.8088)

2 views (last 30 days)
Nbin=? ;
f=? ; %Frequency resolution
EEG_1Spec=? ;
figure;
subplot(2,1,1);
stem(EEG_1Spec);
plot(EEG_1Spec);
title('EEG_1Spec');
xlabel('time(s)');
ylabel('Amplitude(microV)');
subplot(2,1,2);
VEPconvAveSpec=fft(?);
stem(VEPconvAveSpec);
plot(VEPconvAveSpec);
title('VEPconvAveSpec');
xlabel('time(s)');
ylabel('Amplitude(microV)');

Answers (2)

Bora Eryilmaz
Bora Eryilmaz on 22 Dec 2022
Edited: Bora Eryilmaz on 22 Dec 2022
load('VEPdata1.mat')
subplot(211)
plot(EEGdata)
subplot(212)
pspectrum(EEGdata, fs)
ax = gca;
ax.XLim = [0 0.072];
  1 Comment
nazmican
nazmican on 22 Dec 2022
% Before I load VEPdata.1
% I try this code
Nbin=length(t) ;
df=fs/N ; %Frequency resolution
EEG_1Spec=0:70.87;
figure;
subplot(2,1,1);
stem(EEG_1Spec);
plot(EEG_1Spec);
title('EEG_1Spec');
xlabel('frequency(Hz)');
ylabel('Amplitude(microV)');
subplot(2,1,2);
VEPconvAveSpec=fft(127);
stem(VEPconvAveSpec);
plot(VEPconvAveSpec);
title('VEPconvAveSpec');
xlabel('frequency(Hz)');
ylabel('Amplitude(microV)');

Sign in to comment.


Mathieu NOE
Mathieu NOE on 23 Dec 2022
hello
try this
this code split the entire signal is smaller chuncks of NFFT length, do the fft and average the spectrum
also if you are interested only up to 70 Hz , then you can downsample (decimate) your signal from 2000 to 200 Hz
the choice of NFFT depends of what is most important between a finer frequency resolution (higher NFFT) or more averages to reduce noise effects (lower NFFT) ;here I think the optimal value for NFFT lies between 1000 and 2000
Now the peaks in the spectrum appear more distinctively as we have a good frequency resolution with a fair amount of averages
% freq display range
flow = 0;
fhigh = 71;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
filename = 'VEPdata1.mat';
data = load(filename);
signal = (data.EEGdata)';
Fs = data.fs; % sampling rate is TBD Hz
dt = 1/Fs;
[samples,channels] = size(signal);
% time vector
time = (0:samples-1)*dt;
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 10;
if decim>1
for ck = 1:channels
newsignal(:,ck) = decimate(signal(:,ck),decim);
Fs = Fs/decim;
end
signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)*1/Fs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1000; % for maximal resolution take NFFT = signal length (but then no averaging possible)
OVERLAP = 0.75; % overlap will be used only if NFFT is shorter than signal length
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),
plot(time,signal,'b');grid on
title(['Time plot / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
figure(2),plot(freq,20*log10(spectrum_raw),'b');grid on
df = freq(2)-freq(1); % frequency resolution
title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');
xlim([flow fhigh]);
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 percentage of overlap % (between 0 and 0.95)
[samples,channels] = size(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
s_tmp = zeros(nfft,channels);
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*ones(1,channels));
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

Categories

Find more on EEG/MEG/ECoG in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!