Borrar filtros
Borrar filtros

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

3 visualizaciones (últimos 30 días)
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)');

Respuestas (2)

Mathieu NOE
Mathieu NOE el 23 de Dic. de 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

Bora Eryilmaz
Bora Eryilmaz el 22 de Dic. de 2022
Editada: Bora Eryilmaz el 22 de Dic. de 2022
load('VEPdata1.mat')
subplot(211)
plot(EEGdata)
subplot(212)
pspectrum(EEGdata, fs)
ax = gca;
ax.XLim = [0 0.072];
  1 comentario
nazmican
nazmican el 22 de Dic. de 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)');

Iniciar sesión para comentar.

Categorías

Más información sobre EEG/MEG/ECoG en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by