Computation of Power of a Signal in Matlab – Simulation and Verification

1 Star2 Stars3 Stars4 Stars5 Stars (4 votes, average: 5.00 out of 5)

Calculating the energy and power of a signal was discussed in one of the previous posts. Here, we will verify the method by using Discrete Fourier Transform (DFT) in Matlab. Check here to know more on the concept of power and energy.

The total power of a signal can be computed using the following equation
$$P_x=\lim_{N\rightarrow \infty }\frac{1}{N}\sum_{n=0}^{n=N-1}\left | x(n) \right |^2 \rightarrow (1) $$

For other forms of equations: refer here.

Case Study:

\(x(t)\) is a sine wave of amplitude \(A\) and frequency \(f_c\) represented by the following equation.
$$ x(t) = A sin(2 \pi f_c t) $$
When represented in frequency domain, it will look like the one on the right side plot in the following figure. This is evident from the fact that the sinewave can be mathematically represented by applying Euler’s formula.
$$ Asin(2 \pi f_c t) = A\frac{e^{j 2 \pi f_c t}-e^{-j 2 \pi f_c t}}{2j} $$

Taking the Fourier transform of \(x(t)\) to represent it in frequency domain,
$$X(f)=F \left[Asin(2 \pi f_c t) \right]=\int_{-\infty }^{\infty }\left[\frac{e^{j 2 \pi f_c t}-e^{-j 2 \pi f_c t}}{2j}\right]e^{-j2 \pi f t} dt = \frac{A}{2j}\left[ \delta (f-f_c) – \delta (f+f_c) \right ]$$

When considering the amplitude part,the above decomposition gives two spikes of amplitude A/2 on either side of the frequency domain at \(f_c\) and \(-f_c\).

Time domain representation and frequency spectrum of a pure sinewave
Squaring the amplitudes gives the magnitude of power of the individual spikes/frequency components. The power spectrum is plotted below.

Power spectrum of a pure sinewave
Thus if the pure sinewave is of amplitude A=1V and frequency=100Hz, the power spectrum will have two spikes of value \(A^2/4=0.25 \) at 100 Hz and -100 Hz frequencies. The total power will be \(A^2/4+A^2/4=0.25+0.25=0.5 W \)

Let’s verify this through simulation.

Simulation and Verification

A sine wave of 100 Hz frequency and amplitude 1V is taken for the experiment.

A sinusoidal wave of 10 cycles is plotted here

Sine wave generated in Matlab

Matlab’s Norm function:

Matlab’s basic installation comes with “norm” function. The p-norm in Matlab is computed as
$$NORM(v,p) = \left( \sum_{n=0}^{n=N-1}\left | x(v) \right |^p \right)^{1/p} \rightarrow (2)$$

By default, the single argument norm function computed 2-norm given as
$$NORM(v) =NORM(v,2)= \left( \sum_{n=0}^{n=N-1}\left | x(v) \right |^2 \right)^{1/2} \rightarrow (3)$$

To compute the total power of the signal \(x[n]\) (as in equation (1) above), all we have to do is – compute norm(x), square it and divide by the length of the signal.

The above given piece of code will result in the following output

Verifying the total Power by DFT : Frequency Domain

Here, the total power is verified by applying DFT on the sinusoidal sequence. The sinusoidal sequence \(x[n]\) is represented in frequency domain \(X[f]\) using Matlab’s FFT function. The power associated with each frequency point is computed as

$$P_x[f]= X[f]X^{*}[f]$$

Finally, the total power is calculated as the sum of all the points in the frequency domain representation.


Power spectrum of a pure sinewave simulated in Matlab
A word on Matlab’s FFT: Matlab’s FFT is optimized for faster performance if the transform length is a power of 2. The following snippet of code simply calls “fft” without the transform length. In this case, the window length and the transform length are the same. This is to simplify the calculation of power. You can re-write the call to the FFT routine with transform length set to next power of two which is greater than or equal to the window length (sequence length). Then the step to compute total power will be differing slightly in the denominator. But that will not improve resolution (Remember : zero padding to compute FFT will not improve resolution).

Also note that in the above simulation we are using a pure sinusoid. The entire sequence of sinusoid defined all the cycles completely. There are no discontinuities in the sequence. If you call FFT with the transform length set to next power of 2 (as given in Matlab manuals), it will pad additional zeros to the sequence and creates discontinuities in the FFT computation. This will lead to spectral leakage. FFT and spectral leakage is discussed here.

Recommended Books on Signal Processing

  • Mahdi Zaman

    Hello Sear;
    My work about mitigate of nonlinearity on PM-16QAM and PM-QPSK. I dont have any programs about my work like biterate error, OSNR, Spectrum efficiency …etc
    plz if you have any information to help for.
    thanks alot

  • Amey Naik

    To get a quick view of frequency response of a signal we can use,
    freqz(sig,1,length(sig),’whole’) %% gives frequency response on scale of 0 to 2. same as fft normalized on scale of 0 to 2.

  • yasser

    i generate 64 QAM pkt then pass it thru awgn with snr param set to 10 db then generate another 16 QAM pkt pass it thru awgn with same snr of 10 db the problem is when applying the norm method i get different values for the case of 16 and 64 QAM so why is that although they passwd thru same snr of 10 db

    • This method calculates the power of a given signal. It cannot differentiate between signal and noise.
      SNR is a ratio. For your case, the average power of generated 64-QAM may be different from 16-QAM points. Even though they are passed through same SNR, proportional amount of noise will be generated to make SNR=10 dB.

      Since SNR is a ratio, Signal power and noise power can be different for 64/16 QAM packets but maintaining the same SNR.
      Thus adding signal and noise power for both the case will give different result.

    • yasser

      if (ch==0)


      sig=randi([0 1],no_samples,1);

      coded_data =convolution(sig,code_rate,10); %not built in fn in matlab code rate 3/4

      [ x ,scale ]=Modulation(coded_data,mt);%not built in fn in matlab
      pkt64=x*scale ; %for modulation normalization
      info64=awgn(pkt64 , snr , ‘measured’)

      if (ch==1)

      sig=randi([0 1],no_samples,1);

      coded_data =convolution(sig,code_rate,10); %not built in fn in matlab code rate 3/4

      [ x,scale ]=Modulation(coded_data,mt);%not built in fn in matlab
      pkt16=x*scale;%for modulation normalization
      info16=awgn(pkt16 , snr , ‘measured’);

      u see here that symbols of both cases pkt16 , pkt64 are normalized to unity

      • Not sure what’s the Modulation function is returning in place for scale. Normalizing the amplitude to unity ?
        Also the no_samples will be different in both the cases as it depends on M. Try it with a long data sequence. A short sequence will give skewed results.
        Check both the FFT and norm methods

  • yasser

    please i need a recommended method by which in matlab i can estimate SNR of received signal blindly irrespective to modulation scheme /order used

    • Some estimators for AWGN channel can be found here

      Pauluzzi, D.R.; Beaulieu, N.C., “A comparison of SNR estimation techniques for the AWGN channel,” Communications, IEEE Transactions on , vol.48, no.10, pp.1681,1691, Oct 2000
      doi: 10.1109/26.871393
      keywords: {AWGN channels;mean square error methods;noise;parameter estimation;phase shift keying;receivers;AWGN channel;CRB;Cramer-Rao bound;PSK;SNR estimation techniques;baseband 8-PSK signals;binary phase-shift keying;complex AWGN.;estimation techniques;mean square error;performance;real additive white Gaussian noise;signal-to noise ratio;AWGN channels;Additive white noise;Baseband;Computer simulation;Gaussian noise;Mean square error methods;Phase estimation;Phase shift keying;Signal processing;Signal to noise ratio},

      • yasser

        the paper investigates constant modulus modulation tech as MPSK but what about MQAM (non constant modulus) ?

  • yasser

    Px=X.*conj(X)/(L^2) what is the need for L^2

    • It is used for properly scaling FFT output when using L-point FFT/IFFT in Matlab
      This has to be scaled by a factor of 1/L, the length of the sequence.

      Thus for computing power, which is X.*conj(X) (multiplying twice), the scaling factor is 1/(L^2)

      • yasser

        i thought the correct scaling for preserving parseval rule is (length(fft) /sqrt (length(used fft bins = length(data) ) )

        • Several scaling terms are possible.
          1. Scale FFT by dt and IFFT by Fs
          2. Scale FFT by 1/M and IFFT by M
          3. Scale FFT by 1/sqrt(M) and IFFT by sqrt(M)
          4. Scale FFT by 1 and IFFT by 1

          Only important thing here is that the scaling terms in FFT and IFFT will multiply to unity. It also depends on what scaling factor goes on in the implementation.
          You can test with various settings and check which factor gives you the expected result.