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

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)');

Answer 1

[Mathieu NOE]

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

[Mathieu NOE]

hello

Problem solved ?

Answer 2

[Bora Eryilmaz]

load('VEPdata1.mat')

subplot(211)
plot(EEGdata)

subplot(212)
pspectrum(EEGdata, fs)
ax = gca;
ax.XLim = [0 0.072];

[nazmican]

% 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)');