Compute signal power in Matlab

Calculating the energy and power of a signal was discussed in one of the previous posts. Here, we will verify the calculation of signal power 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 = \displaystyle{ \lim_{ N \rightarrow \infty} \sum_{ n=0 }^{ N-1 } | x(n) |^2} \quad\quad (1)

For other forms of equations: refer here.

This article is part of the following books
Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885
Digital Modulations using Python ISBN: 978-1712321638
Wireless communication systems in Matlab ISBN: 979-8648350779
All books available in ebook (PDF) and Paperback formats

Case Study:

Let x(t) be a sine wave of amplitude A and frequency fc represented by the following equation.

x(t) = A\; sin \left( 2 \pi f_c t\right) \quad\quad (2)

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.

A\; sin \left( 2 \pi f_c t\right) = \displaystyle{ A \frac{e^{j 2 \pi f_c t} - e^{-j 2 \pi f_c t} }{2 j}} \quad \quad (3)

Taking the Fourier transform of x(t) to represent it in frequency domain,

\begin{aligned} X(f) &= F \left[ A\; sin \left( 2 \pi f_c t \right)\right] \\ &= \int_{- \infty}^{\infty} \left[  \frac{e^{j 2 \pi f_c t} - e^{-j 2 \pi f_c t} }{2 j}  \right] e^{- j 2 \pi f t} \\ &= \frac{A}{ 2 j } \left[ \delta (f - f_c) - \delta (f + f_c) \right] \end{aligned} \quad\quad (4)

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

Time domain representation and frequency spectrum of a pure sinewave
Figure 1: 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
Figure 2: 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=1; %Amplitude of sine wave
Fc=100; %Frequency of sine wave
Fs=1000; %Sampling rate - oversampled by the rate of 10
Ts=1/Fs; %Sampling period
nCycles=200; %Number of cycles of the sinewave

subplot(2,1,1);
t=0:Ts:nCycles/Fc-Ts; %Time base
x=A*sin(2*pi*Fc*t); %Sinusoidal function
stem(t,x); %Plot command

A sinusoidal wave of 10 cycles is plotted here

Sine wave generated in Matlab
Figure 3: 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( \displaystyle{ \sum_{n=0}^{N-1} |x (v)|^p } \right) ^{\frac{1}{p}} \quad\quad (5)

By default, the single argument norm function computed 2-norm given as

NORM (v) = \left( \displaystyle{ \sum_{n=0}^{N-1} |x (v)|^2 } \right) ^{\frac{1}{2}} \quad\quad (6)

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.

L=length(x);
P=(norm(x)^2)/L;
sprintf('Power of the Signal from Time domain %f',P);

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

Power of the Signal from Time domain 0.500000

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^\ast [f]

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

X = fft(x);
Px=sum(X.*conj(X))/(L^2); %Compute power with proper scaling.
subplot(2,1,2)
% Plot single-sided amplitude spectrum.
stem(Px);
sprintf('Total Power of the Signal from DFT %f',P);

Result:

Total Power of the Signal from DFT 0.500000
Power spectrum of a pure sinewave simulated in Matlab
Figure 4: 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.

Rate this article: PoorBelow averageAverageGoodExcellent (16 votes, average: 4.06 out of 5)

References:

[1] Sanjay Lall,”Norm and Vector spaces”,Information Systems Laboratory,Stanford University.↗

Books by the author

Wireless Communication Systems in Matlab
Wireless Communication Systems in Matlab
Second Edition(PDF)

PoorBelow averageAverageGoodExcellent (162 votes, average: 3.78 out of 5)

Digital modulations using Python
Digital Modulations using Python
(PDF ebook)

PoorBelow averageAverageGoodExcellent (123 votes, average: 3.60 out of 5)

digital_modulations_using_matlab_book_cover
Digital Modulations using Matlab
(PDF ebook)

PoorBelow averageAverageGoodExcellent (126 votes, average: 3.70 out of 5)

Hand-picked Best books on Communication Engineering
Best books on Signal Processing

15 thoughts on “Compute signal power in Matlab”

  1. 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

    Reply
  2. 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.

    Reply
  3. 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

    Reply
    • 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.

      Reply
    • if (ch==0)

      snr=10;
      M=64;
      mt=’64QAM’;
      code_rate=1/2;
      no_samples=.75*fft_size*log2(M)*code_rate;

      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’)
      end

      if (ch==1)

      snr=10;
      M=16;
      mt=’16QAM’;
      code_rate=1/2;
      no_samples=.75*fft_size*log2(M)*code_rate;
      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’);
      end

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

      Reply
      • 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

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

    Reply
    • 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},
      URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=871393&isnumber=18877

      Reply
    • It is used for properly scaling FFT output when using L-point FFT/IFFT in Matlab
      X=fft(x);
      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)

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

        Reply
        • 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
          etc..,

          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.

          Reply

Post your valuable comments !!!