从功率谱生成热图，Matlab

Question

我有一个非常大的局部场电位（原始电压）数据集，我已经对其进行了预处理以消除噪声和异常值。我排列了数据，使每一行代表 30 秒的样本。我生成的功率谱如下：

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

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

目标：生成热图，使 pxx 的每一行都是生成的热图上的一列，所以我应该在 x 轴上有 1075 个 bin，我希望 Y 轴的频率限制在 0 到 120 Hz 之间。我已经尝试使用 imagesc，但遇到了困难，谢谢。

Answer 1

要绘制结果，您需要做一些事情：

1. 选择与高达 120Hz 的频率对应的列；
2. 转置pxx 以使pxx 的行显示为生成图像的列；
3. 如果您希望使用imagesc 将最高频率显示在顶部，请使用flipud 将数据倒置；
4. 可选择转换为对数分贝刻度；
5. 使用imagesc 或pcolor 绘图；
6. 使用caxis 选择要显示的值的动态范围，这样您就可以在颜色图中获得适当的值分布；
7. 为热图样式选择颜色图，例如 colormap(hot)。

这可以通过：

% 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);

或者，如果您想控制显示的轴：

% 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)");

作为一个例子，你会得到一个类似于

的图表

对于一个简单的缓慢频率变化的音调产生的：

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