Hi there, i just want to design and apply a 2s hanning window with 50% overlapp on my EEG signal in order to segment it , i've tried with @hann command but i couldn't get it

my signal is a 17-channel signal with 30721 samples totally (=17 * 30721) . i need my signal to be segmented with 2s hanning window with 50% overlapp to apply preproccessing on it for my Master thesis. sampling frequency is 256 Hz.

i have problem from A to Z designing the window and also, i dont know how to apply it on my Original Signal,

I will be grateful if you could help me design and apply this window to my signal,

thank you so much

Answer 1

[Mathieu NOE]

hello

see example below

the signal gets splitted in 2s buffer and windowed

see the subfunction at the bottom of the code :

        sw = signal(start:stop).*window; % signal buffered and windowed

then we do fft on each windowed segment and when iteration process is complete we get a time / frequency (spectrgram) plot

code :

clearvars

% dummy signal
Fs = 256;
dt = 1/Fs;
t = 0:dt:30;                    % 2 s at 1 kHz sample rate
signal = chirp(t,0,15,50,'quadratic');             % Start at DC, cross 50 Hz at 15 s 
 
nfft = 2*256; % frame length = 2 seconds x Fs
overlap = 0.5; % window overlap ;  0 = 0 % , 0.5 = 50% , 1 = 100 % 

[S,F,T] = myspecgram(signal, Fs, nfft, overlap); 

figure(1); 
imagesc(T,F,20*log10(abs(S))) 
colorbar('vert');
set(gca,'YDir','Normal') 
xlabel('Time (secs)') 
ylabel('Freq (Hz)') 
title('Short-time Fourier Transform spectrum (dB Scale)')
colormap('jet');

function  [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)
% FFT peak spectrogram of signal  (example sinus amplitude 1   = 0 dB after fft).
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    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_specgram = [];
    for ci=1:spectnum
        start = 1+(ci-1)*offset;
        stop = start+nfft-1;
        sw = signal(start:stop).*window; % signal buffered and windowed
        fft_specgram = [fft_specgram abs(fft(sw))*4/nfft];     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end

% 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_specgram = fft_specgram(select,:);
freq_vector = (select - 1)*Fs/nfft;

% time vector 
% time stamps are defined in the middle of the buffer
time = ((0:spectnum-1)*offset + round(nfft/2))/Fs;

end

[alireza ghavami]

Hi dear My friend,

Im so happy to hear from you,

i just saw your answer and immediately, i tried on in my MATLAB , but unfortunately, since im very new to Matlab i could not be successful.

i copied your Function text and the rest in to their own text box in matlab but the result is not the one i need. by the way, i would thank you for your help.

i have an EEG signal and i have to apply preprocess(such as 50 Hz Notch filter in order to denoising electricity noise and bandpass filter[Fc1=0.4 Hz Fc2=45 Hz] to remove other artifacts and i need to apply window for denoising as i think,

this is my problem!! designing and applying Hanning window with 2 second length and 50 % overlapp(which means 512 samples length and 256 samples overlapp).

i've designed the filters but i do not know how to apply them on my signal whether after windowing or without it ...

[Mathieu NOE]

hello again

have you started some code ? I can help you in putting it corectly

do not forget the data as well

tx

[alireza ghavami]

yes dear , i have started to write ...

i have a .mat file which contains 17 channels with 30721 samples each(= 17 * 30721 ). i seperated each channel for preprocessing. i will write my full code but I'm not sure how much of the code is correct...

this is my code:

clc

clear all

close all

load('C:\Users\Alireza\Desktop\DATA MATLAB\data_motad\AbasAshrafi_Close.mat')

%%%%%%%%%%%%%%%%%%%%% Seperating Channels %%%%%%%%%%%%%%%%%%%%

y1 = Close(1,:);

y2 = Close(2,:);

y3 = Close(3,:);

y4 = Close(4,:);

y5 = Close(5,:);

y6 = Close(6,:);

y7 = Close(7,:);

y8 = Close(8,:);

y9 = Close(9,:);

y10 = Close(10,:);

y11 = Close(11,:);

y12 = Close(12,:);

y13 = Close(13,:);

y14 = Close(14,:);

y15 = Close(15,:);

y16 = Close(16,:);

y17 = Close(17,:);

%%%%%%%%%%%%%%%%% Normalizing Signals %%%%%%%%%%%%%%%%

y1_Norm = y1/max(y1);

y2_Norm = y2/max(y2);

y3_Norm = y3/max(y3);

y4_Norm = y4/max(y4);

y5_Norm = y5/max(y5);

y6_Norm = y6/max(y6);

y7_Norm = y7/max(y7);

y8_Norm = y8/max(y8);

y9_Norm = y9/max(y9);

y10_Norm = y10/max(y10);

y11_Norm = y11/max(y11);

y12_Norm = y12/max(y12);

y13_Norm = y13/max(y13);

y14_Norm = y14/max(y14);

y15_Norm = y15/max(y15);

y16_Norm = y16/max(y16);

y17_Norm = y17/max(y17);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fs = 256; %% Sampling Frequency

Nsig = 17; %% number of EEG signal channels

EndTime = size(Close,2)/fs; %% time duration of EEG signal in Sec

Nsamp = EndTime * fs; %% number of samples of EEG signal

times = 1/fs : 1/fs : EndTime; %% time points in sec

%%%%%%%%%%%%%% Segmentation for Applying Filters %%%%%%%%%%%%%%

Segment1 = y1_Norm .* frame;

Segment2 = y2_Norm .* frame;

Segment3 = y3_Norm .* frame;

Segment4 = y4_Norm .* frame;

Segment5 = y5_Norm .* frame;

Segment6 = y6_Norm .* frame;

Segment7 = y7_Norm .* frame;

Segment8 = y8_Norm .* frame;

Segment9 = y9_Norm .* frame;

Segment10 = y10_Norm .* frame;

Segment11 = y11_Norm .* frame;

Segment12 = y12_Norm .* frame;

Segment13 = y13_Norm .* frame;

Segment14 = y14_Norm .* frame;

Segment15 = y15_Norm .* frame;

Segment16 = y16_Norm .* frame;

Segment17 = y17_Norm .* frame;

t = [0:1:30720];

%%%%%%%%%%%%%% Desigining a 50 Hertz Notch Filter %%%%%%%%%%%%%%%%%%%%%

CT_notch_freq = 50;

DT_notch_freq = 2*pi*CT_notch_freq/fs;

% in the Z-plane , zeros = exp(±j*DT_notch_freq) , ...

% and poles = r * exp(±j*DT_notch_freq)

% r is the radius of IIR poles , 0<r<1

r = 0.7;

Notch_zeros=[exp(1i*DT_notch_freq) exp(-1i*DT_notch_freq)];

Notch_poles=[r*exp(1i*DT_notch_freq) r*exp(-1i*DT_notch_freq)];

b_Notch=poly(Notch_zeros);

a_Notch=poly(Notch_poles);

[db , mag , pha , grd , w] = freqz_m(b_Notch,a_Notch);

figure; plot(w * fs / (2 * pi) , db);

%%%%%% Designing a Butterworth Bandpass Filter with [0.4 45] cutoff freq.

n = 6;

Wn = [0.4 45]/(fs/2);

[b_Bandpass , a_BandPass] = butter(n,Wn);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Mathieu NOE]

hello again

can you share your data file as well ? attach your mat file with the paperclip button

[alireza ghavami]

sure, i will attach my file here in this comment

by the way , thank you so much dear.

i've attached my EEG data here, but since it contains 17 channels of data, i think it must be seperated first as i did it on my own code. its name is AbasAshrafi_Close.mat and since its obvious, the data has been recorded when the person had closed his eyes.

[Mathieu NOE]

hello again

so this is a code that will :

1. load your data (time + 16 channels of data) ,
2. do the bandpass and notch filter
3. then split in chuncks of 2s buffers ,
4. apply a window,
5. and finally perform a fft on those buffers and do a linear averaving of the fft outputs

now I am surprised by the look of your data , I am not really a EEG expert but I didn't expect this type of damped sinusoids

when you do the fft spectra of those signals , you can see below how the notch filter removes the 50 Hz line

code :

clc
clearvars

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

filename = 'AbasAshrafi_Close.mat'; 
data = load(filename);
data = (data.Close)'; % let's put the data in column oriented direction
time = data(:,1);  % time vector 
signal = data(:,2:end); % 16 channels of data

dt = mean(diff(time)); % sampling time interval
Fs = 1/dt; % sampling rate is 256 Hz here

% select which channel(s) to process
channels_selected = 8:16;  % can be simply just one value (8) or an array of multiple channels like (8:16) will display data for channels 8 to 16

signal = signal(:,channels_selected);
[samples,channels] = size(signal);

% create some legend char arrays
for ck =1:numel(channels_selected)
    str_leg{ck} = ['Channel :' num2str(channels_selected(ck))];
end

% time vector (this will make sure the time vector has truly constant
% increment (sometimes rounding occurs in excel formated data)
time = (0:samples-1)*dt;

%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
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;

%% band pass  filter section %%%%%%
%bandpass filter[Fc1=0.4 Hz   Fc2=45 Hz]
f_low = 0.4;
f_high = 45;
N_bpf = 2;           % filter order
[b,a] = butter(N_bpf,2/Fs*[f_low f_high]);
signal_filtered = filtfilt(b,a,signal);

%% notch filter 
fc_notch = 50;
w0 = 2*pi*fc_notch/Fs;
p = 0.95; % something slightly less than 1 (will change the width and the depth of the notch filter).

% digital notch (IIR)
num1z=[1 -2*cos(w0) 1];
den1z=[1 -2*p*cos(w0) p^2];

% now let's filter the signal 
signal_filtered = filtfilt(num1z,den1z,signal_filtered);
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 256*2;    % 2 s of data => 2 x Fs samples
OVERLAP = 0.5; % (50% ) overlap will be used only if NFFT is shorter than signal length

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),
subplot(2,1,1),plot(time,signal);grid on
title(['Raw data Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend(str_leg)
subplot(2,1,2),plot(time,signal_filtered);grid on
title(['Filtered data Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend(str_leg)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
[freq, spectrum_filtered] = myfft_peak(signal_filtered,Fs,NFFT,OVERLAP);
df = freq(2)-freq(1); % frequency resolution 

figure(2),
subplot(2,1,1),plot(freq,20*log10(spectrum_raw));grid on
title(['Raw data Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');
legend(str_leg)
subplot(2,1,2),plot(freq,20*log10(spectrum_filtered));grid on
title(['Filtered data Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');
legend(str_leg)


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

[alireza ghavami]

Thank you dear friend...

really i dont know how to appreciate your help

this code is all i need i think, which is the most generous help from you.

if the admin of the website allows me to get in touch with you directly in private, would the best thing for me.

i have to try running my data using your code and then i will tell you the result

thaaaaank youuuu soo much my friend

by the way, Im from Iran and any kind of appreciation i should do to answer a part of your big help, do not hesitate.

[Mathieu NOE]

hello again

my pleasure !

do you mind accepting my answer ?

tx

do not hesitate to come back to me for further assistance

[alireza ghavami]

ohh so sorry .. i accepted now so sorry

thank you so much dear Mathieu

i dont know how to appreciate....

[Mathieu NOE]

no problem

don't be sorry , there are so many times people don't even care to say thank you even when you have spent hours on their problems.

[alireza ghavami]

thank you so much dear Mathieu,

to be honest, i never forget your help. it helped me alot.

I've sent the result to my Supervisor and she said its in time domain but we need frequency domain to see whether it has been denoised or not.

is she right? the result was in time or frequency domain?

if so, how can i fix it?

just one small help, will be the best thing you can do for me. really i dont want to bother you Mathieu, im so rude sorry

[Mathieu NOE]

hello again

see above in my previous answer , you see the two plots ?

the first one is the time domain display (damped sinusoids) and the second one is the frequency domain (fft spectrum)

see the plots x axis label (the first plot x axis is time , on the second plot it's frequency)

I also emphasized that the notch filter has clearly removed the 50 Hz electrical perturbation

Glad I could help you

[alireza ghavami]

Hello again dear Mathieu really i dont know how to appreciate you...you helped me alot and i dont know how to thank you.

my supervisor has agreed me and she said the plots are correct but she wants to see the plot of power spectrum in frequency domain.

during this time, i've searched alot in youtube and other websites to see whether i can find something or not(i'm really ashamed to ask you again and again) but i could not find anything suitable for it.

could you please help me for the last time dear?

[Mathieu NOE]

hello

the power spectrum plot is what I already showed you above

here again

it's coming from this fft computation in my code

[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);

but maybe your supervisor would like you to use this matlab function instead

Analyze signals in the frequency and time-frequency domains - MATLAB pspectrum - MathWorks France

look at the doc and example it's fairly easy to replace my fft code with this matlab function

[alireza ghavami]

hi dear Mathieu , im so sorry to reply you late,

im very thankful for your help, you helped me alot

i was trying to read about your last comment in the help of MATLAB but i did not get it clearly,

and if its possible, can i ask you to write me your contact details so i could contact you easier?

im so sorry if im asking too much questions :)