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

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

[RJDB_BHT]

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));

Answer 1

[Jeremy]

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);

[Sudhar Ekanayake]

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;

Answer 2

[Sudhar Ekanayake]

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;