使用 Python 或 Matlab，我如何创建一个直方图，其中每个 bin 等于一个成比例的频率范围（例如一个八度音阶），并且添加到每个 bin 的数据是在该频率范围内具有任何幅度的样本数量的计数?
这是我到目前为止所拥有的：
length = data.shape[0] / Fs
print(f"length = {length}s")
RMS = lambda x: np.sqrt(np.mean(x**2))
sample = np.arange(int(length))
RMS_of_sample = np.zeros(sample.shape)
for ns in sample:
RMS_of_sample[ns] = RMS(data[ns*Fs:(ns+1)*Fs])
plt.hist(RMS_of_sample, label="Left channel")
plt.show()
但这似乎创建了一个等于整个光谱范围的 bin，因此直方图中只出现了一个条形。我也不确定条形图中的数据是否真的代表幅度出现的计数。
如果您想要确定信号在特定频带中的能量大小，您可以对信号执行 FFT。为此，您需要确保您有振幅数据而不是能量数据。所以你会想要声压Pa而不是Pa^2。在计算 FFT 之后，这些值必须乘以才能得到 RMS 格式。
根据您的时间数据长度和采样时间，您可以获得频域中信号的特定频率分辨率。使用个时间样本和的采样时间，您将获得 Herz 的频率分辨率。现在，如果您想知道特定频率范围的能量贡献，您可以使用之前以 RMS 格式计算的线性频谱来计算该频率范围的 RMS 值。根据 Parseval 定理，可以使用以下等式计算特定频率范围的 RMS 值：
其中 s 是频谱中不同频率的幅度：
我使用这种方法来计算一个信号的 10 Hz 频段的能量，该信号包含两个幅度为且频率和 Hz 的正弦波。采样频率为 0.0098。
这是时间信号：
产生的具有峰峰值的线性光谱为：
从中我计算了 Hz 箱的 RMS 值。所以第一条表示 Hz 信号的 RMS 值。第二个从 Hz 依此类推。
如果添加了一些噪声，则光谱变为：
Hz 分档的 RMS 值变为：
我进一步建议查看均方根和总体水平
这是我用来执行上述操作的 Matlab 代码（尽管可能仍然包含错误）：
clear all
close all
t=linspace(0,10,1024); % time vector
y=2*sin(5*2*pi*t)+2*sin(45*2*pi*t)+randn(1,length(t)); % time signal
Fs=1/(t(2)-t(1));
Fn=Fs/2;
N=length(y);
yf=fft(y,N);
df = Fs/N; % Frequency Resolution
plot(t,y)
grid on
title('Time Signal')
ylabel('Amplitude')
xlabel(['Time Resolution: ',num2str(1/Fs),' s'])
amplH = abs(yf);
amplitudengang = fftshift(amplH/N);
x_fn = 0 : df : Fn-df;
x_fa = 0 : df : Fs-df;
figure
stem(x_fa-Fn, amplitudengang, 'b.-')
axis([-Fn Fn 0 max(amplitudengang)])
title('Two sided Spectra')
ylabel('Amplitude')
xlabel(['Frequency Resolution: ',num2str(df),' Hz'])
grid
amplitudengang=[amplitudengang(513) amplitudengang(514:end).*2];
figure
stem([0:df:(Fn-df)], amplitudengang, 'b.-')
axis([0 Fn 0 max(amplitudengang)])
title('Single Sided Spectra')
ylabel('Amplitude')
xlabel(['Frequency Resolution: ',num2str(df),'Hz'])
grid
amp10=amplitudengang(1:round((10-df)/df));
amp20=amplitudengang(round((10-df)/df):round((20-df)/df));
amp30=amplitudengang(round((20-df)/df):round((30-df)/df));
amp40=amplitudengang(round((30-df)/df):round((40-df)/df));
amp50=amplitudengang(round((40-df)/df):round((50-df)/df));
rms10=amp10./sqrt(2);
rms20=amp20./sqrt(2);
rms30=amp30./sqrt(2);
rms40=amp40./sqrt(2);
rms50=amp50./sqrt(2);
rms10=sqrt((rms10(1)/2)^2+sum(rms10(2:(end-1)).^2)+(rms10(end)/2)^2);
rms20=sqrt((rms20(1)/2)^2+sum(rms20(2:(end-1)).^2)+(rms20(end)/2)^2);
rms30=sqrt((rms30(1)/2)^2+sum(rms30(2:(end-1)).^2)+(rms30(end)/2)^2);
rms40=sqrt((rms40(1)/2)^2+sum(rms40(2:(end-1)).^2)+(rms40(end)/2)^2);
rms50=sqrt((rms50(1)/2)^2+sum(rms50(2:(end-1)).^2)+(rms50(end)/2)^2);
figure
stem([5 15 25 35 45],[rms10,rms20,rms30,rms40,rms50])
axis([0 Fn 0 max([rms10,rms20,rms30,rms40,rms50])])
title('RMS values per 10Hz')
ylabel('Amplitude')
xlabel(['Frequency in Hz'])
grid on
编辑
为了减少泄漏，需要对时间数据进行窗口化。之后，生成的光谱需要乘以匹配的校正因子：
有两种不同的方法可以做到这一点：
除非您有大量的波段，否则方法 1 是更好的方法。在线查找一些倍频程滤波器或自己设计。每个频段都有一个滤波器。要获得 1 kHz 频段的能量，只需使用 1 kHz 倍频程滤波器对输入进行滤波。
然后将“运行”RMS 应用于输出。有很多不同的方法可以调整时域行为。最简单的方法是将其分割成帧并计算每帧的 RMS。您还可以做更多花哨的事情：非常快速的启动和缓慢的衰减，这是典型的“米桥”会做的。
对于 FFT 方法，您可以将信号切成重叠的帧，应用一个窗口，然后进行 FFT。然后为每个八度音程（或第三个八度音程等）应用不同的频谱掩码，并计算每个频带中的（复数）RMS。光谱掩模应通过查看 ANSI S1.11-2004 得出，这是“可接受”滤光片的标准。这种方法在调整输出的时域属性方面的灵活性较低，您只需为每个波段和每帧获得一个 RMS 值。