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

조회 수: 150 (최근 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 CenterFile 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