FFT vs DFT Issue with Power Calculation

I am trying to compare MATLAB fft with Mike’s DFT Complex Sine Wave Signal dot product. Because DFT loop is very slow it takes too long to perform practical calculations. So I use a scaling factor. But when I use this factor, the size of my power spectrum results are smaller than fft results. Likewise when my scale factor is 1, the comparison between fft and DFT seem exactly the same.

Am I doing something wrong in calculating power spectrum from DFT?

%%
%     PURPOSE: Compare fft with dft for an audio file
%  AUDIO FILE: Use a public domain birdcall from xeno-canto.org
%      RESULT: Produce 3 plots: 1) Time domain, 2) FFT, and 3) DFT
%   DEVELOPER: http://www.biophysicslab.com/portfolio/matlab-projects/
%        DATE: June 6, 2020
%
%  REFERENCES:
%   "Signal Processing Problems on Udemy.com" 
%       by Dr. Mike X Cohen, sigprocMXC_SpectBirdcall.m
%   "Master the Fourier transform and its applications on Udemy.com" 
%       by Dr. Mike X Cohen, Fourier_DTFT.m
%
%%

% Load in birdcall (source: https://www.xeno-canto.org/403881).
[bc,fs] = audioread('XC403881.mp3');

% Configuration params for Fourier transforms.
freq_range = [0 8000]; % Hz
dft_loop_reduction = 100; % set to 1 for no reduction (Hint: takes a long time)

% let's hear it!
% soundsc(bc,fs)

n = length(bc);
hz = linspace(0,fs/2,floor(n/2)+1);
signal = detrend(bc(:,1)); % smooth bird call audio audio for fft & dft

% Create a time vector based on the data sampling rate.
timevec = (0:n-1)/fs;

% Plot the data from the two audio file channels.
figure(1), clf
subplot(311)
% Include a small offset for left and right channels.
plot(timevec,bsxfun(@plus,bc,[.2 0]))
xlabel('Time (sec.)')
title('Time domain')
set(gca,'ytick',[],'xlim',timevec([1 end]))

%% Compute & plot the power spectrum using MAYTLAB fft function.

bcpow_fft = abs(fft( signal )/n).^2;

subplot(312)
plot(hz,bcpow_fft(1:length(hz)),'linew',2)
xlabel('Frequency (Hz)')
ylabel('Power');
title(['Frequency domain using FFT with ' num2str(length(bcpow_fft)) ' points']);
set(gca,'xlim',freq_range)
% Make fft and dft y-limits be the same for comparison.
ylim_fft = get(gca,'ylim'); 

%% Compute & plot the power spectrum using DFT loop.

fourTime = (0:n-1)/n;
fCoefs   = zeros(size(signal));

h = waitbar(0,'Please wait for DFT loop...');
for fi=1:dft_loop_reduction:n
    
    % Create complex sine wave.
    csw = exp( -1i*2*pi*(fi-1)*fourTime );
    
    % Compute dot product between sine wave and signal (Fourier coefficients).
    fCoefs(fi) = sum( signal'.*csw ) / n;
    
    % GUI to show progress for long calculation times.
    waitbar(fi/n)
end
close(h)

bcpow_dft = abs(fCoefs).^2;

subplot(313)
plot(hz,bcpow_dft(1:length(hz)),'linew',2)
xlabel('Frequency (Hz)')
ylabel('Power');
title(['Frequency domain using DFT loop with ' ...
    num2str(floor(length(bcpow_dft)/dft_loop_reduction)) ' points']);
set(gca,'xlim',freq_range)
set(gca,'ylim',ylim_fft)

%% done.

Hey, how many times did I tell you never to use the loop-based DFT??!!

I think this is a visual effect. You’re missing a lot of frequencies in the DFT. Try plotting them as
plot(bcpow_dft,bcpow_fft,'o')
You’ll see a perfect correlation and then lots of zeros on the left. I suspect if you run a correlation on all the downsampled points (that is, the non-zero elements of bcpow_dft), you’ll find the correlation to be >.99 (unlikely to be exact because of precision errors).

btw, when running code in a loop that takes a long time, try to remove anything that doesn’t need to be in the loop. In your case, that includes the transpose, the /n, and the waitbar updating. It’s like factoring out constants from a summation.

1 Like

Thank you Mike.
After taking three of our MATLAB classes on Udemy I thought I would try and answer a question on MATLAB Central for the first time. The question was how to avoid using FFT to get DFT of an audio file. Sadly, just after answering the question, the person who asked the question clobbered his question. But the MATLAB Central “overlord” tried to repair the damage so at least my answer has some context. Never seen this behavior before…

https://www.mathworks.com/matlabcentral/answers/544352-test1-code1-test2-code2#answer_449087

Ah, the internet is a strange place… but kudos for posting there! Indeed, your reply is still visible even if the original question isn’t, so perhaps someone will find it useful.