相位谱的1/3倍频程在matlab中如何画呢？

Question

相位谱的1/3倍频程在matlab中如何画呢？
如下图所示，是一个音频文件的相位谱和幅值谱，其中图1是没有倍频程，图2是有倍频程，如何实现相位的倍频程呢？

Answer 1

平头小混混 上午好☀️☀️☀️️
本答案参考通义千问

要绘制Matlab中的相位谱及其1/3倍频程，您可以使用以下步骤：

1. 首先，计算原始相位谱：

   % 假设您已经有了音频信号 'audioSignal' 和采样率 'fs'
   [Pxx,Fs] = pwelch(audioSignal,[],[],[],Fs); % 计算幅度谱 (Pxx) 和采样率 (Fs)
   phs = unwrap(angle(Pxx)); % 提取相位数据（unwrap函数用于处理相位跳跃）
   

2. 然后，计算1/3倍频程边界频率：

   f半个octave = Fs / 3; % 每个半音阶对应的频率间隔
   freqBins = linspace(0, Fs/2, length(phs)/2+1); % 相位谱频率范围
   octaveBands = floor(freqBins / f半个octave) * f半个octave; % 半音阶带宽的频率点
   
   % 创建一个与相位谱相同长度的向量来存储每个1/3倍频段的起始和结束相位
   phaseBoundaries = zeros(size(octaveBands));
   for i = 1:length(octaveBands)-1
       phaseBoundaries(i*2-1) = min(find(freqBins >= octaveBands(i)));
       phaseBoundaries(i*2) = max(find(freqBins <= octaveBands(i+1)));
   end
   

3. 绘制原始相位谱和1/3倍频程分界线：

   figure;
   subplot(2,1,1);
   plot(freqBins, phs, '-r', 'LineWidth', 1.5); % 原始相位谱
   hold on;
   plot([freqBins(phaseBoundaries), freqBins(phaseBoundaries)], [-180, 180], '--g'); % 1/3倍频程分界线
   legend('Phase Spectrum', '1/3 Octave Boundaries');
   xlabel('Frequency (Hz)');
   ylabel('Phase (deg)');
   title('Original Phase Spectrum');
   
   subplot(2,1,2);
   % 绘制带有1/3倍频程分界的相位谱，可以根据需要选择显示方式
   surf(freqBins, ones(length(freqBins))/length(freqBins), phs);
   caxis([-180, 180]);
   view(-90, 90); % 设置视角以查看整个频谱
   hold on;
   patch(freqBins(phaseBoundaries), ones(size(phaseBoundaries))/length(freqBins), phaseBoundaries, 'FaceColor', 'g', 'EdgeColor', 'none');
   legend('Phase Spectrum', '1/3 Octave Boundaries');
   xlabel('Frequency (Hz)');
   ylabel('Bin Index');
   zlabel('Phase (deg)');
   title('Phase Spectrum with 1/3 Octave Bands');
   

确保将 'audioSignal' 替换为您的实际音频信号，并设置正确的采样率 Fs。执行以上代码后，您将看到一个包含原始相位谱及其1/3倍频程分界线的图像。