MATLAB Answers

Plotting spectrogram for EEG data

113 views (last 30 days)
Sparsha Kumari
Sparsha Kumari on 5 Feb 2020
Commented: Walter Roberson on 8 Feb 2021 at 1:36
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.

  3 Comments

Jakob B. Nielsen
Jakob B. Nielsen on 5 Feb 2020
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
Sparsha Kumari on 5 Feb 2020
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
Sudon Biswas on 6 Feb 2021 at 9:54
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

Sign in to comment.

Answers (2)

Daniel M
Daniel M on 24 Apr 2020
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.

  5 Comments

Show 2 older comments
Sudon Biswas
Sudon Biswas on 6 Feb 2021 at 21:06
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
Sudon Biswas on 6 Feb 2021 at 21:07
@Walter Roberson
Walter Roberson
Walter Roberson on 7 Feb 2021 at 4:27
Coherence is only defined between two signals, not one signal by itself. https://en.wikipedia.org/wiki/Coherence_(signal_processing)

Sign in to comment.


Sudon Biswas
Sudon Biswas on 6 Feb 2021 at 21:22
Edited: Sudon Biswas on 6 Feb 2021 at 21:23
@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

  3 Comments

Walter Roberson
Walter Roberson on 7 Feb 2021 at 4:29
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
Babu Biswas on 7 Feb 2021 at 16:09
%% %% 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
Walter Roberson on 8 Feb 2021 at 1:36
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.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by