RMS power in signal: time series vs frequency series. PSD/RMS definition?

109 次查看（过去 30 天）
RJDB_BHT
RJDB_BHT 2013-7-16
评论： Sudhar Ekanayake 2018-7-8
I am building a function that computes the PSD and RMS power in a time series. Idon't have a license to the DSP library.
Below my function, and a test file. When ran 'as is', it computes the RMS power based on the time series (variable: avPow) and the RMS power using the PSD (variable: r). For a simple (complete) sine wave, it return the correct value of 0.5*sqrt(2).
With the other test signal it can be shown that the amplitudes in the fft are correct.
My question: apparently, with PSD = (signalFFT(1:n) .* signalFFT(1:n)); RMS = 0.5 * sqrt(sum(PSD)/n*Fs)
variable r returns the correct result. Why is there a factor 0.5?
Any help is much appreciated.
Test file:
clear clc close format compact
Fs = 1000; % sample frequency t = 0:1/Fs:1; % time x = sin(2*pi*t); % > rms of sin(x) = 1/sqrt(2)
% x = 1 + sin(2.5*pi*t) + 2*cos(0.5*pi*t); % x = 2*cos(2*pi*t*10)+cos(2*pi*t*20); % x = 0.3*randn(size(t, 1), size(t, 2)); % plot(t, x)
n = size(t, 2); % number of data points avPow = sqrt(1 / n * sum((x - mean(x)) .* (x - mean(x)))); % average power based on time series [f, p, r] = computePSD(x, 1/Fs, false); % f: frequencies; p: power; r: rms power
disp(avPow) disp(r)
Function: function [ F, PSD, RMS ] = computePSD( signal, sampleTime, plotFigs ) [r, c] = size(signal); L = max(r, c);
if (c ~= 1 && r ~= 1)
error('encountered non-row or non-column array');
end
% sample frequency = 1 / sampleTime
Fs = 1 / sampleTime;
% a signal with time interval dF can contain frequency components with a
% frequency of at most Fs / 2 -> Nyquist frequency
F = Fs * (0:(L/2)) / L;
[r2, c2] = size(F);
n = max( r2, c2);
signal = signal - mean(signal);
signalFFT = 2 * abs(fft(signal)) / L;
PSD = (signalFFT(1:n) .* signalFFT(1:n));
% not a RMS in the regular sense; the area under the PSD is the square of
% the RMS noise, as according to http:%www.youtube.com/watch?v=ywChrIRIXWQ
%RMS = sqrt(sum((abs(signalFFT).^2)/length(n)/Fs));
RMS = 0.5 * sqrt(sum(PSD)/n*Fs);
if (plotFigs)
% Plot single-sided amplitude spectrum.
figure
loglog(F, signalFFT(1:n));
% dummy = gca();
% dummy.log_flags = 'lln';
title('FFT of signal');
xlabel('Frequency (Hz)');
ylabel('|signal(f)|');
figure
loglog(F, PSD)
% dummy = gca();
%dummy.log_flags = 'lln';
title('PSD of signal');
xlabel('Frequency (Hz)')
ylabel('|signal^2 / Hz|')
end
end
  1 个评论
RJDB_BHT
RJDB_BHT 2013-7-16
编辑：RJDB_BHT 2013-7-16
Additionally, I noticed that when I change the time from 0:1 to 0:tmax in the test files, r underestimates the RMS power by a factor sqrt(tmax)
This yields the same answer for the power computed with the time series and the frequency series: RMS = sqrt(0.5*sum(PSD));

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Jeremy
Jeremy 2015-6-18
Your equations are a little muddled, and there are three basic errors.
  1. There is a factor of 2 on the your fft result, I assume this is for the conversion from a double sided spectra to a single sided spectra. This factor should be on the PSD and NOT the fft and so this increases the PSD by a factor of 2 (rms by a factor of 1.44)
  2. Your PSD clalculation should also be divided by the bin width and this is why your result is changing when you use a longer time series (bin width is no longer 1.0)
  3. In the RMS calulation you multiply the sum of the PSD by Fs/n (2.0). the RMS is sqrt of the integral of the PSD so it should be the sum of the PSD times the bin width which is 1.0. This error accounts for the other factor of 2 on the PSD (1.44 for the RMS). see the corrected lines below from your function:
signal = signal - mean(signal);
signalFFT = abs(fft(signal)) / L;
binWidth=Fs/length(signal); %almost exactly 1 with 1.001 sec time history.
PSD = 2*(signalFFT(1:n) .* signalFFT(1:n))/(binWidth);
% not a RMS in the regular sense; the area under the PSD is the square of
% the RMS noise, as according to http:%www.youtube.com/watch?v=ywChrIRIXWQ
%RMS = sqrt(sum((abs(signalFFT).^2)/length(n)/Fs));
RMS =sqrt(sum(PSD)*binWidth);
  1 个评论
Sudhar Ekanayake
Sudhar Ekanayake 2018-7-8
RMS =sqrt(sum(PSD)*binWidth); suppose to divide by signal length, L for the mean value ? or is it done at PSD calculation signalFFT = 2 * abs(fft(signal)) / L;

请先登录，再进行评论。

Sudhar Ekanayake
Sudhar Ekanayake 2018-7-1
RMS =sqrt(sum(PSD)*binWidth); suppose to divide by signal length, L for the mean value ? or is it done at PSD calculation signalFFT = 2 * abs(fft(signal)) / L;

类别

在 Help Center 和 File Exchange 中查找有关 Parametric Spectral Estimation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by