pwelch函数计算功率出错

% 参数设置

Fs = 1e6; % 采样率 1 MHz

t_end = 1; % 信号时长 1秒

t = 0:1/Fs:t_end-1/Fs; % 时间向量

f1 = 1e3; % 1 kHz 信号

f2 = 2e3; % 2 kHz 信号

f3 = 5e3; % 5 kHz 信号

f4 = 13e3; % 13 kHz 信号

% 生成信号

signal1 = 1*sin(2*pi*f1*t); % 1 kHz 正弦波，幅值1

signal2 = 1*sin(2*pi*f2*t); % 2 kHz 正弦波，幅值1

signal3 = 1*sin(2*pi*f3*t); % 5 kHz 正弦波，幅值1

signal4 = 1*sin(2*pi*f4*t); % 13 kHz 正弦波，幅值1

% 合成信号

combined_signal = signal1 + signal2 + signal3 + signal4;

% PSD功率谱计算

[psd, f] = pwelch(combined_signal, [], [], 4096*4096, Fs); % 计算功率谱密度

f_range = f >= 20 & f <= 20000; % 限制频率范围为 20-20000 Hz

psd_range = psd(f_range);

f_range = f(f_range);

plot(f_range, psd_range);

xlabel('Frequency (Hz)');

ylabel('Power Spectral Density (dB/Hz)');

title('Power Spectral Density');

grid on;

sum(abs(combined_signal).^2)/length(combined_signal) %信号功率

P=Fs/2/length(psd_range)*(sum(psd_range)-0.5*(psd_range(1)+ psd_range(end)))%梯形法积分

ans =

2

P =

50.0500

Answer 1

[UDAYA PEDDIRAJU]

Hi ZW,

I did slight modification to your code, it is working fine for me. Converting "psd_range" to decibels (dB) for better visualization.

Let me know if there is anything specific you want to know?

% Parameter settings
Fs = 1e6; % Sampling rate 1 MHz
t_end = 1; % Signal duration 1 second
t = 0:1/Fs:t_end-1/Fs; % Time vector

% Frequencies
frequencies = [1e3, 2e3, 5e3, 13e3]; % Frequencies of the sine waves
amplitude = 1; % Amplitude of the sine waves

% Generate composite signal
combined_signal = sum(amplitude * sin(2 * pi * frequencies' * t), 1);

% PSD power spectrum calculation
[psd, f] = pwelch(combined_signal, [], [], 4096*4096, Fs); % Calculate power spectrum density

% Limit the frequency range to 20-20000 Hz
f_range = f >= 20 & f <= 20000; 
psd_range = psd(f_range);
f_range = f(f_range);

% Plotting
figure;
plot(f_range, 10*log10(psd_range)); % Convert to dB
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB/Hz)');
title('Power Spectral Density');
grid on;

% Calculate signal power
signal_power = sum(abs(combined_signal).^2) / length(combined_signal);
P = Fs / 2 / length(psd_range) * (sum(psd_range) - 0.5 * (psd_range(1) + psd_range(end)));