Plotting spectrogram for EEG data

Greetings!

I am done with the filtering and plotting the power spectrum for the EEG data. Now I want to plot the spectrogram for this data. I have tried using the 'spectrogram' and 'stft' function and played with the various input options but I do not get the plot the way I want it to be. Could someone please show me some light in this regard.

Thanks!

P.S. : I need to get it as shown in figure 1 (Reference: William E. DeCoteau et. al. www.pnas.orgcgidoi10.1073pnas.0700818104). Instead I get it as shown in the figure 2.

[Jakob B. Nielsen]

Can you show us what you do to get your figure, and your data? Otherwise, off the top of my mind, try using the surf function and add shading interp in there.

[Sparsha Kumari]

Thanks a lot for your reply! I get it as shown in figure 2. I have used the 'stft' function as follows and got the plot,

stft(data, fs);

[Sudon Biswas]

Could you plz provide me spectrogram and magnitude squared coherence (MSC) Matlab code for ECG signal? The code must be for single signal example: 1. For raw signal 2. Noisy Signal 3. For noise free signal

Answer 1

[Daniel M]

You should plot the absolute value of your output from spectrogram or stft. Maybe even the log, or 10*log10 too (but I don't recall exactly, you'll have to look up the documentation). Basically, when you call those functions with outputs, try and make the plot match the one you get when you call the functions without outputs.

Your spectral and temporal resolution are very different from the first figure. That paper is doing a time-frequency analysis over ~ 4 seconds and 25 Hz. You are looking at ~150 ms and 500 Hz. You can look at this tutorial (with video) on how time window and frequency resolution trade-off with each other.

http://www.fieldtriptoolbox.org/workshop/oslo2019/timefrequency/

[Sudon Biswas]

Could you plz provide me spectrogram and magnitude squared coherence (MSC) Matlab code for ECG signal? The code must be for single signal example: 1. For raw signal 2. Noisy Signal 3. For noise free signal

[Walter Roberson]

@Sudon Biswas

Your wording sounds a lot like you have a homework assignment and you are hoping that someone will write the code for you instead of you writing the code yourself.

[Sudon Biswas]

Actually I am trying but can not understand the clearly thats why I am asking Because in matlab command there are two input signal.but I want to try for each signal MSC Sir

[Sudon Biswas]

@Walter Roberson

[Walter Roberson]

Coherence is only defined between two signals, not one signal by itself. https://en.wikipedia.org/wiki/Coherence_(signal_processing)

Answer 2

[Sudon Biswas]

@Walter Roberson

Sir, Here is the code

clc; load('100m.mat'); %100.m data from physionet ECG = (val-0)/200; fs = 360; %sampling frequency t = (0:length(ECG)-1)/fs; %time period(total sample/Fs ) plot(t,ECG, 'k'); xlabel('Samples'); ylabel(' Amplitude [mV] '); title( 'Heart beat' ); hold on; % legend('ECG Signal'); % %%figure(); %Time domain ploting % subplot(411); plot(t,ECG, 'k'); axis([0 10 -1.5 1.5]); grid on; xlabel(' Samples '); ylabel(' Amplitude [mV] '); title(' Original ECG Signal '); %Time domain ploting %figure();

1. * Item one
2. * Item two

%% %% Generationg Power line interference (noise) signal = ECG; %original ecg signal fNoise = 50; % Frequency [Hz] aNoise = 0.20; % Amplitude noise = aNoise*sin(2*pi.*t.*fNoise); % subplot(412); plot(t, noise, 'k'); % %legend('PLI'); grid on; axis([0 10 -1.5 1.5]); xlabel(' Time [s] '); ylabel(' Amplitude [mV] '); % title(' 50Hz Power Line Intereference (PLI) '); %%figure();

%% %% 50Hz Noise adding with the Original ECG Signal Noisy_Signal = (signal + noise); % subplot(412); plot(t, Noisy_Signal, 'k'); grid on; axis([0 10 -1.5 1.5]); xlabel(' Samples '); ylabel(' Amplitude [mV] '); title(' ECG Signal with 50Hz PLI'); % figure();

%% %%Spectral analysis of the ECG signal with 50Hz noise (PLI)signal L = length(Noisy_Signal); NFFT = 2^nextpow2(L); noise_fft = abs(fft(Noisy_Signal,NFFT)); %%% creating frequency axis freq = fs/2*linspace(0,1,NFFT/2+1);

%%%%% plot single-sided amplitude spectrum of % subplot(2,3,4); subplot(411); plot(freq,noise_fft(1:NFFT/2+1)); grid on; axis([0 500 0 500]); title('Single-Sided Amplitude Spectrum of ECG Signal with 50Hz PLI'); xlabel('Frequency (Hz)'); ylabel('|Amplitude, Y(f)|'); % figure();

% % %% % % FIR Low_pass filter design fs1=1000; fc=100; Wc = fc/(fs1/2); %cut_off frequency N = 40; % N= order of filter b = fir1(N,Wc,'low',rectwin(N+1)); % Rectangular window % % freqz(b,1,1024,fs1); % frequency response of FIR LPF filter in Hz % % freqz(b,1,512); % frequency response of FIR LPF filter in radian. % % fvtool(b,1); % % %%figure(); % %% %% filter the noisy signal f_fir = filter(b,1,Noisy_Signal); % subplot(513); % plot(t,f_fir); % reconstructed_Signal = f_fir; % grid on; % title(' Rectangular window filtered signal '); % Rectangular window % xlabel(' Samples '); % ylabel(' Amplitude [mV]'); % axis([0 10 -1.5 1.5]); %%figure()

Actually I am trying but can not understand the clearly thats why I am asking Because in matlab command there are two input signal.but I want to try for each signal MSC 1. Original ECG signal 2. Noisy_ECG signal 3. Filtered ECG Signal

[Walter Roberson]

Please replace the above with a formatted copy of your code. Go into the editor and delete the code you posted. Then click the ">" button in the header to get a code entry box. Go over to MATLAB and copy your code. Return to your browser and Paste the code into the code entry block.

[Babu Biswas]

%% %%    ECG signal generating from MIT-BIH arithmia database
clear all;
close all;
clc; 
load('100m.mat');           %100.m data from physionet
ECG = (val-0)/200; 
Fs = 360;                    %sampling frequency
t = (0:length(ECG)-1)/Fs;    %time period(total sample/Fs )  
plot(t,ECG);    
xlabel('Samples');
ylabel(' Amplitude [mV] ');
title( 'Heart beat' );
hold on;
% legend('ECG Signal');
figure();                                      %Time domain ploting
% subplot(411); 

plot(t,ECG);
axis([0 10 -1.5 1.5]);
grid on;
xlabel(' Samples ');
ylabel(' Amplitude [mV] ');
title(' Original ECG Signal ');                 %Time domain ploting
figure();
%%  %%  Generationg Power line interference (noise) 
signal = ECG;     %original ecg signal
fNoise = 50;    % Frequency [Hz]
aNoise = 0.20;  % Amplitude 
noise  = aNoise*sin(2*pi.*t.*fNoise);
% subplot(412); 
plot(t, noise);
% %legend('PLI');
grid on;
axis([0 10 -1.5 1.5]);
xlabel(' Time [s] ');
ylabel(' Amplitude [mV] ');
% title(' 50Hz Power Line Intereference (PLI) ');
figure();

%%  %% 50Hz Noise adding with the Original ECG Signal
Noisy_Signal = (signal + noise);
% spectrogram(Noisy_Signal,'yaxis')
pspectrum(Noisy_Signal,Fs,'spectrogram','TimeResolution',0.5)
title('Normal Signal')
% subplot(412); 
plot(t, Noisy_Signal);
grid on;
axis([0 10 -1.5 1.5]);
xlabel(' Samples ');
ylabel(' Amplitude [mV] ');
title(' ECG Signal with 50Hz PLI');
% figure();
%%   %%Spectral analysis of the ECG signal with 50Hz noise (PLI)signal
L = length(Noisy_Signal);
NFFT = 2^nextpow2(L);
noise_fft = abs(fft(Noisy_Signal,NFFT));
%%%  creating frequency axis
freq = Fs/2*linspace(0,1,NFFT/2+1);

%%%%% plot single-sided amplitude spectrum of 
% subplot(2,3,4);
% subplot(411);
plot(freq,noise_fft(1:NFFT/2+1)); 
grid on;
axis([0 500 0 500]);
title('Single-Sided Amplitude Spectrum of ECG Signal with 50Hz PLI');
xlabel('Frequency (Hz)');
ylabel('|Amplitude, Y(f)|');
% figure();

% % %%  % % FIR Low_pass  filter design
fs1=1000;
fc=100;
Wc = fc/(fs1/2);                               %cut_off frequency
N = 40;                                        % N= order of filter
b = fir1(N,Wc,'low',rectwin(N+1));             % Rectangular window
% % freqz(b,1,1024,fs1);                       % frequency response of FIR LPF filter in Hz
% % freqz(b,1,512);                            % frequency response of FIR LPF filter in radian.
% % fvtool(b,1);
% figure();
% 
%% %% filter the noisy signal
f_fir = filter(b,1,Noisy_Signal);
% subplot(513); 
plot(t,f_fir); 
reconstructed_Signal = f_fir;
grid on;
title(' Rectangular window filtered signal ');          % Rectangular  window
xlabel(' Samples ');
ylabel(' Amplitude [mV]');
axis([0 10 -1.5 1.5]);
% figure()
% 
%%%%% plot single-sided amplitude spectrum of Rectangular window
%%   %%Spectral analysis of the ECG signal with 50Hz noise (PLI)signal
L = length(f_fir);
NFFT = 2^nextpow2(L);
noise_fft = abs(fft(f_fir,NFFT));
%%%  creating frequency axis
freq = Fs/2*linspace(0,1,NFFT/2+1);
% subplot(2,3,4); 
% subplot(412);
plot(freq,noise_fft(1:NFFT/2+1));                  % Rectangular window
grid on;
axis([0 500 0 500]);
title('Rectangular filtered ECG Signal');
xlabel('Frequency (Hz)');
ylabel('|Amplitude, Y(f)|');
figure();

want to know

Actually I am trying but can not understand the clearly thats why I am asking Because in matlab command there are two input signal.but I want to try for each signal MSC

1. Original ECG signal

2. Noisy_ECG signal

3. Filtered ECG Signal

Thanks for your interest for helping me.

[Walter Roberson]

That's fine, but MSC is simply not defined for one signal by itself.

MSC would be defined for comparing the filtered ECG signal to the original ECG signal.