% EE3220 - Dr. Durant - 2014-02-18 % DFT and scaling of power % Summary: Simply taking 20*log10(abs(DFT_value)) gives you relative % decibels that are useful for COMPARING power spectra of 2 signals of the % same length. This approach ignores the reference level. What are the dB % relative to? The code below shows one standard reference level: 0 dB % represents 1 W assuming that the signal voltage is applied across a 1 Ohm % resistive load. There are many other standard reference levels in use; % see http://en.wikipedia.org/wiki/Decibel#Voltage fs = 10; % Sampling frequency, Hz f = 1; % Signal frequency, Hz A = 5; % Amplitude, Volts T = [1 100]; % Duration, s; 2 values to demonstrate the difference between energy and power N = T*fs; % Number of samples for each T t = (0:max(N)-1)/fs; % Time values in s - we'll use a subset when the selected T is not the maximum T x = A * cos(2*pi*f*t); % Voltage waveform figure for idxT = 1:length(T) % Create a graph pair for each duration Y = fft(x(1:N(idxT))); normFreqs = (0:N(idxT)-1) / N(idxT) * 2; % Number of pi's at frequency sampling point subplot(2,length(T),idxT) % 1st plot will be UNNORMALIZED dB plot(normFreqs,20*log10(abs(Y))) title(sprintf('Spectra for %g s sinusoid', T(idxT))) if idxT==1, ylabel('Magnitude, unnormalized (dB)'), end % only label leftmost graph subplot(2,length(T),idxT+length(T)) % 2nd plot will have dB NORMALIZED relative to 1 W assuming 1 Ohm load Y = Y / sqrt(N(idxT)); % Correction #1: DFT increases energy; this is corrected by 1/N factor in IDFT. % To avoid energy scaling, DFT and IDFT each need to be scaled by 1/sqrt(N), % but that requires an expensive square root and divisions. Note that % the 1/N in the IDFT is often just a shift since N is often a power of 2. Y = Y / sqrt(N(idxT)); % Correction #2: DFT now conserves energy, but it is standard to plot % a POWER spectrum, not an ENERGY spectrum. So, we need to divide % energy by time; on the voltage side of the equation, this is % equivalent to dividing by the square root of time. plot(normFreqs,20*log10(abs(Y))) xlabel('Normalized digital frequency (radians/sample \times\pi)') if idxT==1, ylabel('Magnitude, normalized (dB)'), end % only label leftmost graph end % So, was the normalization correct? Let's see Vrms = A/sqrt(2) % RMS voltage given sinusoid amplitude Pavg = Vrms^2 % Average power, W, assuming 1 Ohm load PavgHalf = Pavg/2 % Power is divided between + and - frequencies in DFT dbExpected = 10*log10(PavgHalf) % convert to dB re 1 W