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
For other forms of equations: refer here.
This article is part of the following books |
Case Study:
Let x(t) be a sine wave of amplitude A and frequency fc represented by the following equation.
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.
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)](https://s0.wp.com/latex.php?latex=%5Cbegin%7Baligned%7D+X%28f%29+%26%3D+F+%5Cleft%5B+A%5C%3B+sin+%5Cleft%28+2+%5Cpi+f_c+t+%5Cright%29%5Cright%5D+%5C%5C+%26%3D+%5Cint_%7B-+%5Cinfty%7D%5E%7B%5Cinfty%7D+%5Cleft%5B++%5Cfrac%7Be%5E%7Bj+2+%5Cpi+f_c+t%7D+-+e%5E%7B-j+2+%5Cpi+f_c+t%7D+%7D%7B2+j%7D++%5Cright%5D+e%5E%7B-+j+2+%5Cpi+f+t%7D+%5C%5C+%26%3D+%5Cfrac%7BA%7D%7B+2+j+%7D+%5Cleft%5B+%5Cdelta+%28f+-+f_c%29+-+%5Cdelta+%28f+%2B+f_c%29+%5Cright%5D+%5Cend%7Baligned%7D+%5Cquad%5Cquad+%284%29++&bg=ffffff&fg=000&s=2&c=20201002)
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.
Squaring the amplitudes gives the magnitude of power of the individual spikes/frequency components. The power spectrum is plotted below.
Thus if the pure sinewave is of amplitude A=1V and frequency=100Hz, the power spectrum will have two spikes of value 
at 100 Hz and -100 Hz frequencies. The total power will be 
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 commandA sinusoidal wave of 10 cycles is plotted here
Matlab’s Norm function:
Matlab’s basic installation comes with “norm” function. The p-norm in Matlab is computed as
By default, the single argument norm function computed 2-norm given as
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.500000Verifying 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]](https://s0.wp.com/latex.php?latex=P_x+%5Bf%5D+%3D+X%5Bf%5D+X%5E%5Cast+%5Bf%5D++&bg=ffffff&fg=000&s=2&c=20201002)
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
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: 
(17 votes, average: 3.88 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 Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. |  Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. |  Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. |
| Hand-picked Best books on Communication Engineering Best books on Signal Processing |
||
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
I don’t have it either. sorry.
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.
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.
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
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
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},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=871393&isnumber=18877
the paper investigates constant modulus modulation tech as MPSK but what about MQAM (non constant modulus) ?
Please search IEEE website
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
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)
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
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.