spectrogram强制使用Hz为频率单位

通常，我们使用matlab内置函数spectrogram来绘制信号的时频图，非常的方便。
但是对于下图，显然我们只关心100Hz以下的频段。

时频图

使用ylim([0 0.1])命令可以控制纵轴频显示的频率范围。但是频率以kHz为单位，看着很不方便。

调整频率显示范围
修改频率轴以Hz为单位的操作如下：

方法1 假修改
假修改是指欺骗眼睛，实际上图中实际坐标并没有更改。如下图Datatip所示，实际Y轴坐标没有改变，只是改变了yticklabel。

yticks(0:0.01:0.1)
yticklabels(num2cell(0:10:100))
ylabel('Frequency (Hz)')

或者另一段代码可以有同样的效果，自动化的程度更高些。

ax = gca;
yRuler = ax.YAxis;
% Loop through TickValues, multiply and insert into TickLabels cell-array
for n = 1:numel(yRuler.TickValues)
    yRuler.TickLabels{n} = num2str(yRuler.TickValues(n) * 1000);
end

假修改

方法2 真修改
真修改是指图中坐标与坐标轴的显示是一致的。

• step 1 双击spectrogram字样，右键打开spectrogram函数；
• step 2 在函数内找到第168行“pspectrogram({x},'spect',varargin{:});”，双击pspectrogram字样，右键打开pspectrogram函数；
• step 3 在函数第163、164行，有两句“[f,~,uf] = engunits(f,'unicode');”、“freqlbl = getfreqlbl([uf 'Hz']);” ，这两句控制着频率显示的单位。将这两句替换为“freqlbl = getfreqlbl([ 'Hz']);”，以后使用spectrogram就可以强制使用Hz为频率单位了。

% [f,~,uf] = engunits(f,'unicode');% 自动调整频率单位
% freqlbl = getfreqlbl([uf 'Hz']);  % 自动调整频率单位
freqlbl = getfreqlbl([ 'Hz']);   % 强制使用Hz，禁止自动调整频率单位

• step 4 建议不直接修改内置函数，重命名一个自己的函数会更好，让原函数保留自动调整频率范围的功能。除非你愿意改来改去。

真修改

下面附上我修改好的绘图函数程序代码以及调用格式。

[S,f,t]=spectrogram(wt,4096,4096-16,4096,fs_w,'yaxis');
figure;displayspectrogram(t,f,S)

function displayspectrogram(t,f,Pxx,faxisloc,esttype,thresh)

% 本函数修改自matlab自带同名函数，用于绘制STFT时频图
% 2019/10/24 Leung
% faxisloc  填写'yaxis'或'xaxis'或者空着，默认yaxis，即横轴时间，纵轴频率
% esttype   填写{'psd','power'}或者空着，默认'psd'
% thresh    sets to zero those elements of ps such that 10 log10(ps) ≤ thresh. 
%           Specify thresh in decibels.
%           即设置colorbar的下界。如-20dB应设为：-20
% see "help spectrogram" for complete description of all input arguments.

narginchk(3,6); % narginchk(minArgs,maxArgs)

switch nargin
    case 3
        faxisloc = 'yaxis';
        esttype = 'psd';
        threshold = 0;
    case 4
        esttype = 'psd';
        threshold = 0;
    case 5
        threshold = 0;
    case 6
        threshold = 10^(thresh/10);
end
        
% remove low-power estimates if requested
if threshold>0
  Pxx(Pxx<threshold) = 0;
end


% Cell array of the standard frequency units strings
% Use engineering units
% [f,~,uf] = engunits(f,'unicode');% 自动调整频率单位
% freqlbl = getfreqlbl([uf 'Hz']);  % 自动调整频率单位
freqlbl = getfreqlbl([ 'Hz']);   % 强制使用Hz，禁止自动调整频率单位
[t,~,ut] = engunits(t,'unicode','time');
timelbl = [getString(message('signal:spectrogram:Time')) ' (' ut ')'];


h = newplot;
if strcmpi(faxisloc,'yaxis')
  xlbl = timelbl;
  ylbl = freqlbl;
else
  xlbl = freqlbl;
  ylbl = timelbl;
end

hRotate = uigettool(ancestor(h,'Figure'),'Exploration.Rotate');
if isempty(hRotate) || strcmp(hRotate.State,'off')
  if strcmp(faxisloc,'yaxis')
    hndl = imagesc(t, f, 10*log10(abs(Pxx)+eps));
  else
    hndl = imagesc(f, t, 10*log10(abs(Pxx)'+eps));
  end
  hndl.Parent.YDir = 'normal';

  setupListeners(hndl);
else
  if strcmp(faxisloc,'yaxis')
    hndl = surf(t, f, 10*log10(abs(Pxx)+eps),'EdgeColor','none');
  else
    hndl = surf(f, t, 10*log10(abs(Pxx)'+eps),'EdgeColor','none');
  end
  axis xy
  axis tight
  view(0,90);
end  

if threshold>0
  Pmax = max(Pxx(:));
  if threshold < Pmax
    set(ancestor(hndl,'axes'),'CLim',10*log10([threshold Pmax]))
  end
end

if strcmpi(esttype,'power')
  cblabel = getString(message('signal:dspdata:dspdata:PowerdB'));
else
  cblabel = getString(message('signal:dspdata:dspdata:PowerfrequencydBHz'));
end
%sigutils.internal.colorbari('titlelong',cblabel);
h = colorbar;
h.Label.String = cblabel;

ylabel(ylbl);
xlabel(xlbl);


% -------------------------------------------------------------------------
function setupListeners(hndl)
hAxes = ancestor(hndl,'Axes');
hRotate = uigettool(ancestor(hndl,'Figure'),'Exploration.Rotate');
eYScale = addlistener(hAxes,'YScale','PreSet',@(src,evt) image2surf(hndl));
eView = addlistener(hAxes,'View','PostSet',@(src,evt) image2surf(hndl));
if ~isempty(hRotate)
    eRotate = addlistener(hRotate,'State','PostSet',@(src,evt) image2surf(hndl));
else
    eRotate = [];
end

if ~isprop(hndl,'TransientUserDataListener')
    pi = addprop(hndl,'TransientUserDataListener');
    pi.Transient = true;
end

set(hndl,'TransientUserDataListener',{eYScale,eView,eRotate});