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

118 views (last 30 days)
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 Comment
RJDB_BHT
RJDB_BHT on 16 Jul 2013
Edited: RJDB_BHT on 16 Jul 2013
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));

Sign in to comment.

Answers (2)

Jeremy
Jeremy on 18 Jun 2015
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 Comment
Sudhar Ekanayake
Sudhar Ekanayake on 8 Jul 2018
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;

Sign in to comment.


Sudhar Ekanayake
Sudhar Ekanayake on 1 Jul 2018
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;

Community Treasure Hunt

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

Start Hunting!