How to create a matrix with power spectral density for a range of frequencies across time

Hello, I currently have a data set from electrophysiological recordings which have voltage over time for different trials. Columns represent time (1ms) and rows represent trials. I would like to convert this matrix to represent power (values) over time (columns) at a range of frequencies (rows). If anyone knows how to do this please let me know!

[Dave B]

@Andy Garcia Barrabeitg - this sounds very similar to your last post:

https://www.mathworks.com/matlabcentral/answers/1569208

Did you run into trouble using spectrogram?

[Andy Garcia Barrabeitg]

Hi Dave I did. The output was not what I am looking for unfortunately. I get a 513x8 complex double. I'd like to retrieve the output dimensions outlined in the above question and then plot it from there. I'm able to retain the dimensions using pwelch but would like to extend this further to find the power at different frequency values. For the previous question, the code I used was as follows. Ideally I'd like to reach an output similar to the attached image

data = NCL0243_349_HC_LFPEventMatrix_A.rawdata; 
numTrials = size(data,1);
% create a 3-D matrix which is [n,m,k] big, n x m are going to be the size 
% of a single spectrogam and k is the number of trials (for us, rows)
for i = 1:numTrials 
      sGram(:,:,i) = spectrogram(data(i,:));
end
%Takes the mean over the third dimension 
meanMatrix = (mean(sGram,3)); %works
plot(meanMatrix)

[Dave B]

I think you had a few issues here:

• The fact that the result was complex
• The fact that the shape didn't 'match up'
• The fact that you got (I'm guessing) a spectrogram with very little time resolution (8 bins, the default of the specrogram function)

You can get the PSD estimate using the 4th arg from spectrogram. Or take abs(s), or do exp(abs(s)) (quite common in neurophys)

It would be weird for the shape to be the same as your input. As I mentioned in the previous post, you can still gather x and y values from the call to spectrogram and pass those in to image.

The spectrogram won't be useful for ephys without passing in at least a window and overlap. You need to think about how big you want your window and how much overlap you want. This depends on how you collected the data and what questions you want to answer.

For example:

win=1000; % How many samples for each fft.
          % More samples means better frequency resolution, less time resolution
          % If your sampling rate was 1kHz, this would mean each spectrum
          % is calculated based on 1 second of data
          
over=900;  % How much overlap between spectrum calculations
           % If your sampling rate was 1kHz, and win is 1000, this means
           % slide the window over by .1s at a time

% A quick way to initialize parameters, rather than computing the size of
% the spectrogram:
[s,f,t]=spectrogram(lfp(1,:), win, over, [], fs); 

% now we have f and t which we'll use later, and s which can be used to
% initialize the matrix that will hold all the spectra
s_all=nan([size(s) 3]);

for i = 1:size(lfp,1)
    s=spectrogram(lfp(i,:), win, over, [], fs);
    s_all(:,:,i)=abs(s); % or use p, if you prefer to abs(s)
    
    % or, if you want psd estimates:
    %[~,~,~,p]=spectrogram(lfp(i,:), win, over, [], fs);
    %p_all(:,:,i) = p;
end

imagesc(t,f,log(mean(s_all,3))) % or log10 if you preft
axis xy
ylim([0 60])
colorbar
colormap jet

[Andy Garcia Barrabeitg]

Hi Dave this is almost perfect, thank you! Just two things, could you show me how to implement the log10? I tried log10(abs(s)) as you'd shown in the previous post but that didn't work, I receive the error "Dimension argument must be a positive integer scalar, a vector of unique positive integers, or 'all'". This may naturally fix the next issue which is the expected power values would be between -55 to -25. Never forced them to be this way with NeuroExplorer just doing a methods check and hoping it aligns with the MATLAB ouptut.

[Dave B]

Hmmm, I'd think it would be log10(mean(s_all,3))? I couldn't decide if it should be the mean of the log or the log of the mean. But either way I don't know why you'd run into an issue with log10(abs(s)) - the error you described looks like something I'd expect if you passed a second argument to mean...?

Once you get it to work, you may find that you need 10*log10 as this would be decibels - not sure why that's an appropriate unit though (or why any unit is particularly appropriate for the types of questions neuroscience normally asks of these sorts of data)

[Andy Garcia Barrabeitg]

Hi Dave,

That helped get that going, thank you. I'm still seeing that the vector f is not directly representative of frequency. I'm hoping to have a final matrix with columns representing time (good so far) and the rows representing frequency with the values representing power (good so far). A separate issue I noticed is that the time seems to start at an offset, which I can minimize by having small windows and big overlaps but the output then starts looking odd and it's ultimately not dealing with the foundational issue.

[Dave B]

The vector f should represent frequency to the extent that you gave spectorgram the correct sampling rate. For the vector t: it's not that it starts at an offset so much as each set of values refers to a window. For example, if you had no overlap and you window size was 1 second, your first spectrum is for the time 0-1s.

Here's a demo of both.

fs=20000;
t=linspace(0,3,fs*3);
y(1:fs)         = sin(t(1:fs)*1000*2*pi);
y(fs+1:2*fs)    = sin(t(1:fs)*2500*2*pi);
y(2*fs+1:3*fs)  = sin(t(1:fs)*7777*2*pi);
[s,f,t]=spectrogram(y, fs, 0, [], fs);
imagesc(t,f,log(abs(s)))
colormap hot
caxis([0 1])
axis xy
yline(1000,'b','LineWidth',1)
yline(2500,'b','LineWidth',1)
yline(7777,'b','LineWidth',1)

[Andy Garcia Barrabeitg]

Hi Dave thanks for that answer that helped a ton. I'm now wondering why the time for the spectrogram is cutting out the first and last half a second of data, do you happen to know why?