如何在 dBFS 中绘制 PSD?

How to plot PSD in dBFS?

我编写的以下代码首先采用 FFT 测量信号 y 的 PSD。

Fsig=26e3;           % signal frequency
Fs = 1e6;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
N = 2^15;             % Length of signal
t = (0:N-1)*T;        % Time vector

y=sin(2*pi*Fsig*t) + harmonics; // a sine function with harmonics

% find FFT+PSD
xdft = fft(y);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/N:Fs/2;
psdx_log=10*log10(psdx);

% plot PSD
figure(1)
plot(frq,psdx_log,'r')
ylim([-150 0])
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

这是输出:

如您所见,除了基频之外,频谱还包含各种谐波。然而,这个情节并不是我想要的。我想以 dBFS 为单位绘制 PSD,其中基音正弦波处于 0 dBFS,以便所有谐波都相对于此 0 dBFS 参考进行测量。怎么办?

您可以对PSD进行归一化,这样PSD的最大值将为0dB(如果某些谐波高于信号,您应该调整该值进行归一化):

Fsig=26e3;           % signal frequency
Fs = 1e6;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
N = 2^15;             % Length of signal
t = (0:N-1)*T;        % Time vector

y=sin(2*pi*Fsig*t) + harmonics; // a sine function with harmonics

% find FFT+PSD
xdft = fft(y);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/N:Fs/2;
psdx_log=10*log10(psdx./max(psdx));

% plot PSD
figure(1)
plot(frq,psdx_log,'r')
ylim([-150 0])
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')