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

3 visualizaciones (últimos 30 días)
nazmican el 22 de Dic. de 2022
Comentada: Mathieu NOE el 31 de En. de 2023
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)');
##### 0 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

### Respuestas (2)

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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
filename = 'VEPdata1.mat';
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
##### 1 comentarioMostrar -1 comentarios más antiguosOcultar -1 comentarios más antiguos
Mathieu NOE el 31 de En. de 2023
hello
Problem solved ?

Iniciar sesión para comentar.

Bora Eryilmaz el 22 de Dic. de 2022
Editada: Bora Eryilmaz el 22 de Dic. de 2022
subplot(211)
plot(EEGdata)
subplot(212)
pspectrum(EEGdata, fs)
ax = gca;
ax.XLim = [0 0.072];
##### 1 comentarioMostrar -1 comentarios más antiguosOcultar -1 comentarios más antiguos
nazmican el 22 de Dic. de 2022
% 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