How can I combine complex spectrogram results as an image for multi-channel EEG signal?

I would like to put my multi-channel EEG spectrogram results upside down to get a single image and ad save it. How can I do this? Thank you.

[Daniel M]

Don't know what you mean by upside down, but if it's complex, you could just take the abs() of your signal to get the magnitude.

[Funda KUTLU ONAY]

I want to put each spectrogram result under the other one into the image form. For example, if i have 8-channel, i must obtain 8 different spectrogram. I have to get one image by putting each one under another.

[Daniel M]

help subplot

[Funda KUTLU ONAY]

That's exactly what I want. To obtain a single image from different channels.

Or how can I use the spectrogram() function to get such an image?

Answer 1

[Daniel M]

Take all your outputs from spectrogram, and assuming they are all sampled at the same time/freqs, just stack them (concatenate) vertically. Then you can use either imagesc or surf (with view([0,90])) to plot them. They both have their own different characteristics. Then just set your ticks and labels to be how you want them. See this post for how to do that.

https://www.mathworks.com/matlabcentral/answers/487050-how-to-plot-like-the-inserted-picture#answer_398084

[Funda KUTLU ONAY]

Thank you. I have tried, but I can't adjust the colors. The whole image is getting blue. If I send a sample signal to your mail address, would you try?

[Daniel M]

Can you upload a .mat file here, along with the script you're using to create the figure?

[Funda KUTLU ONAY]

Fs=500; ws=250; ov=ws/2; nfft = 2^nextpow2(length(dt));

ss=[];

for i=1:8

[s f t]=spectrogram(dt(:,i),chebwin(ws),ov,Fs);

ss=[ss;s]

end

[s f t]=spectrogram(dt(:,i),chebwin(ws),ov,Fs);

surf(t, f, 20*log10(abs(s)), 'EdgeColor', 'none');

axis xy;

axis tight;

colormap(jet); view(0,90);

Now I can't set the t and f values. So, I can only draw one right now.

How do I draw 'ss' data which is the combination of 8 channels?

Answer 2

[Daniel M]

It usually takes a lot of work to get a figure to look just right. Here is what I came up with quickly, you can play around with it. See the attached figure. I left extra tidbits in the comments that might help your process.

You might want to try smoothing your data, either at the raw data level (i.e. additional preprocessing to remove noise, rereference channels, notch filter, etc.), at the spectrogram level (different windows, different overlaps), or at the visual level (there are functions that can smooth figures like this, such as smoothn on the FileExchange).

clearvars
close all
clc

mydata = load('mydata.mat');
dt = mydata.dt;

% look at the effects of detrending
figure
subplot(2,1,1)
plot(dt)
title('Raw data')
subplot(2,1,2)
plot(detrend(dt))
title('Detrended data')

dt = detrend(mydata.dt); 

Fs=500; 
ws=250; 
ov=ws/2; 
nfft = 2^nextpow2(length(dt));

% I'm sure you could preallocate s based on chebwin(ws) and ov, etc.
for i=1:8
    [s(:,:,i), f, t]=spectrogram(dt(:,i),chebwin(ws),ov,Fs);
    % since f and t should all be equal.
end
% you probably want to specify your own frequency set, as an input to
% spectrogram()

val = 20*log(abs(s));

% figure
% s1 = surf(t,f,val(:,:,1));
% view([0 90])
% set(s1,'EdgeColor','none')
% title('surf')
% xlabel('Time [s]')
% ylabel('Frequency [Hz]')
% c1 = colorbar;
% ylabel(c1,'Power/frequency')
% axis tight

f1 = figure;
i1 = imagesc(t,f,val(:,:,1));
ax1 = gca;
set(ax1,'YDir','normal')
% title('imagesc')
xlabel('Time [s]')
ylabel('Frequency [Hz]')
c2 = colorbar;
ylabel(c2,'Power/frequency')

% concatenate them such that val2 = 
% [val(:,:,1); val(:,:,2);....; val(:,:,end)]
val2 = reshape(permute(val,[1 3 2]),[],19,1);
nchans = size(dt,2);
ff = repmat(f,nchans,1);
ff = 1:size(val2,1);

% get ytick and yticklabel based off ax1
ax1_tick = get(ax1,'YTick');
ftick_ax1 = round(interp1(f,1:length(f),ax1_tick(1:2:end-1))); % rounding is fine here

% need to replicate this 
ftick = ftick_ax1' + length(f).*(0:nchans-1);
ftick = ftick(:);

ax1_label = get(ax1,'YTickLabel');
fticklabel = repmat(ax1_label(1:2:end-1),nchans,1);

% Look at a histogram of val2 to get a sense of the tails. If there are
% some extreme values, your colorbar (showing the min and max) will look
% odd.

% figure
% histogram(val2(:))

% Plot all channels
f2 = figure('Position',[680 378 805 600]);
i2 = imagesc(t,ff,val2);
ax2 = gca;
set(ax2,'YDir','normal')
xlabel('Time [s]')
ylabel('Frequency [Hz]')
c3 = colorbar;
ylabel(c3,'Power/frequency')
set(ax2,'YTick',ftick,'YTickLabel',fticklabel)
Ftsz = get(ax2,'FontSize'); % get original fontsize because yax changes it

yax = get(ax2,'YAxis');
set(yax,'FontSize',12);
set(get(yax,'Label'),'FontSize',Ftsz);

% set caxis based on percentile of val2. 
[prc] = prctile(val2(:),[1,99]);
cmin = prc(1);
cmax = prc(2);

% can also try based on distribution
mu = mean(val2(:));
sd = std(val2(:));
cmin = mu-3*sd;
cmax = mu+3*sd;

caxis([cmin,cmax])

[Arthur Batte]

hello Daniel, thx for this wonderful code. I personally am not a good programmer however, i would like to create something similar to that of Funda but using your code together with the dataset provided by funda seem not to work. it is crashing somewhere at ftick = ftick_ax1'+ length(f).*(0:nchans-1).

Matrix dimensions must agree. how would you go around this?

[Daniel M]

Try ftick_ax1 without the apostrophe

[Arthur Batte]

Even after doing that, the problem is still the same. ftick_ax1 produces 1 81 160 as the output.

length (f).*(0:nchans-1) produces 0 251 502 753 1004 1255 1506 1757 as the output.

as you can see, the dimensions dont match. i cant figure out where the problem is. Kindly advise

[Daniel M]

You must be using an older version of MATLAB, before implicit expansion was the default.

Change the line

ftick = ftick_ax1' + length(f).*(0:nchans-1);

to

ftick = bsxfun(@plus,ftick_ax1',length(f).*(0:nchans-1));

Let me know if there are any other issues.

[Arthur Batte]

ftick_ax1 is undefined. problem is with ftick_ax1 = round(interp1(f,1:length(f),ax1_tick(1:2:end-1)))

[Daniel M]

You had ftick_ax1 before. Did you try it with new data or something? Could you upload it

[Arthur Batte]

I am still using the same code you provided before and the same data provided mydata but for some unknown reason, its not running to the end. it crashes around this line

ftick_ax1 = round(interp1(f,1:length(f),ax1_tick(1:2:end-1))

clearvars close all clc

mydata = load('mydata.mat'); dt = mydata.dt;

% look at the effects of detrending figure subplot(2,1,1) plot(dt) title('Raw data') subplot(2,1,2) plot(detrend(dt)) title('Detrended data')

dt = detrend(mydata.dt);

Fs=500; ws=250; ov=ws/2; nfft = 2^nextpow2(length(dt));

% I'm sure you could preallocate s based on chebwin(ws) and ov, etc. for i=1:8 [s(:,:,i), f, t]=spectrogram(dt(:,i),chebwin(ws),ov,Fs); % since f and t should all be equal. end % you probably want to specify your own frequency set, as an input to % spectrogram()

val = 20*log(abs(s));

% figure % s1 = surf(t,f,val(:,:,1)); % view([0 90]) % set(s1,'EdgeColor','none') % title('surf') % xlabel('Time [s]') % ylabel('Frequency [Hz]') % c1 = colorbar; % ylabel(c1,'Power/frequency') % axis tight

f1 = figure; i1 = imagesc(t,f,val(:,:,1)); ax1 = gca; set(ax1,'YDir','normal') % title('imagesc') xlabel('Time [s]') ylabel('Frequency [Hz]') c2 = colorbar; ylabel(c2,'Power/frequency')

% concatenate them such that val2 = % [val(:,:,1); val(:,:,2);....; val(:,:,end)] val2 = reshape(permute(val,[1 3 2]),[],19,1); nchans = size(dt,2); ff = repmat(f,nchans,1); ff = 1:size(val2,1);

% get ytick and yticklabel based off ax1 ax1_tick = get(ax1,'YTick'); ftick_ax1 = round(interp1(f,1:length(f),ax1_tick(1:2:end-1))); % rounding is fine here

% need to replicate this ftick = ftick_ax1' + length(f).*(0:nchans-1); ftick = ftick(:);

ax1_label = get(ax1,'YTickLabel'); fticklabel = repmat(ax1_label(1:2:end-1),nchans,1);

% Look at a histogram of val2 to get a sense of the tails. If there are % some extreme values, your colorbar (showing the min and max) will look % odd.

% figure % histogram(val2(:))

% Plot all channels f2 = figure('Position',[680 378 805 600]); i2 = imagesc(t,ff,val2); ax2 = gca; set(ax2,'YDir','normal') xlabel('Time [s]') ylabel('Frequency [Hz]') c3 = colorbar; ylabel(c3,'Power/frequency') set(ax2,'YTick',ftick,'YTickLabel',fticklabel) Ftsz = get(ax2,'FontSize'); % get original fontsize because yax changes it

yax = get(ax2,'YAxis'); set(yax,'FontSize',12); set(get(yax,'Label'),'FontSize',Ftsz);

% set caxis based on percentile of val2. [prc] = prctile(val2(:),[1,99]); cmin = prc(1); cmax = prc(2);

% can also try based on distribution mu = mean(val2(:)); sd = std(val2(:)); cmin = mu-3*sd; cmax = mu+3*sd;

caxis([cmin,cmax])

[Daniel M]

But how did you get it to run before? You managed to get past the line

ftick_ax1 = round(interp1(f,1:length(f),ax1_tick(1:2:end-1))); % rounding is fine here

previously. So what is happening now? What is the exact error message? What version of MATLAB are you using? Are you running this in the command window or in a script?

[Arthur Batte]

Hello Daniel,

Now that you are mentioning the version, it seem that the problem was with the version of matlab. I was initially using the version of 2009.

when i changed it to 2018b, it worked well.

Thank you for the efforts you rendered to me.

[Arthur Batte]

hello Daniel, it appears that the spectrograms created by the code above are stacked horizontally. I would wish to stack them vertically ( see figure attached) where by the horizontal axis is replaced by months of the year and the vertical axis is kept as frequency. In this case, each channel represents a certain month in which it was recorded. For the data given above, they are 8 channels that makes it january to august.

what changes can be do to this code?

[Daniel M]

I recommend that you open a new question regarding this.

Answer 3

[Funda KUTLU ONAY]

Thank you very much for your time. I think I can use it this way. But I have one more question for you. Can I save only the spectrogram part as an image?

[Daniel M]

f2 = figure('Position',[680 378 805 600]);
i2 = imagesc(t,ff,val2);
ax2 = gca;
set(ax2,'YDir','normal')
set(get(ax2,'YAxis'),'Visible','off')
set(get(ax2,'XAxis'),'Visible','off')
caxis([cmin,cmax])