Generating heatmap from power spectrum, Matlab

Question

I have a very large data set of local field potentials (raw voltages) that I have pre-processed to remove noise and outliers. I arranged the data so that each row represents 30 seconds of samples. I have generated power-spectrums as follows:

Fs = 1024
LFP = 1075x30720 double
pxx = 1075x4097 double

for k = 1:1075;
    pxx(kk,:) = pwelch(LFP(k,:));
end

Goal: generate heatmap so that each row of the pxx is a column on generated heatmap, so I should have 1075 bins on the x axis and I'd like to have the Y axis limited to frequencies from 0 - 120 Hz. I've tried using imagesc but am having difficulty, thank you.

Answer 1

To plot the result you would need to do a few things:

1. select columns corresponding to frequencies up to 120Hz;
2. transpose pxx to get the rows of pxx to appear as columns of the generated image;
3. flip the data upside down with flipud if you want the highest frequency to appear at the top with imagesc;
4. optionally convert to logarithmic decibels scale;
5. plot using imagesc or pcolor;
6. choose the dynamic range of values to display with caxis so you get a decent spread of the value in the colormap;
7. select a colormap such as colormap(hot) for a heatmap style.

This can be done with:

% 1) Compute maximum frequency index
Fmax = 120; % Hz
M = 1 + round(Fmax/(0.5*Fs/(size(pxx,2)-1)));
%    select displayed section 
pxx_select  = pxx(:,1:M);

% 2) transpose matrix
pxx_reshape = transpose(pxx_select);

% 3) flip data upside down for imagesc
pxx_reshape = flipud(pxx_reshape);

% 4) convert to decibel scale
pxx_dB      = 10*log10(pxx_reshape);

% 5) generate plot
figure(1);
imagesc(pxx_dB);

% 6) choose dynamic range
% assign e.g. 80dB below peak power to the first value in the colormap
%             and the peak power to the last value in the colormap
caxis([max(max(pxx_dB))-80 max(max(pxx_dB))]);

% 7) select colormap
colormap(hot);

Or, if you want to have control over the displayed axis:

% 1) Compute maximum frequency index
Fmax = 120; % Hz
M = 1 + round(Fmax/(0.5*Fs/(size(pxx,2)-1)));
%    select displayed section 
pxx_select  = pxx(:,1:M);

% 2) transpose matrix
pxx_reshape = transpose(pxx_select);

% 3) flipud not needed with pcolor, instead set t & f axis: 
t = (size(LPF,2)/Fs)*[0:size(LPF,1)];
f = [0:M]*Fmax/(M-1);

% 4) convert to decibel scale
pxx_dB      = 10*log10(pxx_reshape);

% 5) generate plot
figure(2);
% Note: extend by one row & column since pcolor does not show the last row/col
P2 = [pxx_dB pxx_dB(:,end); pxx_dB(end,:) pxx_dB(end,end)];
pcolor(t,f,P2); shading flat;

% 6) choose dynamic range
% assign e.g. 80dB below peak power to the first value in the colormap
%             and the peak power to the last value in the colormap
caxis([max(max(pxx_dB))-80 max(max(pxx_dB))]);

% 7) select colormap
colormap(hot);

xlabel("time (s)");
ylabel("frequency (Hz)");

As an illustration, you would get a graph similar to

for a simple slowly frequency varying tone generated with:

T = size(LPF,1)-1;
phi = 0;
n = [0:size(LPF,2)-1];
for k=1:size(LPF,1)
  f = 0.5*(fmin+fmax) + 0.5*(fmax-fmin)*sin(2*pi*k/T);
  LPF(k,:) = sin(2*pi*f*n/Fs + phi);
  phi = mod(phi + 2*pi*f*size(LPF,2)/Fs, 2*pi);
end