% log approx
clear all; clc
figure(1); clf

% some random data to use for power spectrum
t = 0:0.1:50*pi;
v = sin(t) + 5*randn(size(t));

%plot the waveform
subplot(3,1,1)
 plot(t,v)
 title('time waveform')
 % compare exact and approx FFT mag
 % see also http://en.wikipedia.org/wiki/Alpha_max_plus_beta_min_algorithm
 V = fft(v);
 subplot(3,1,2);
 %exact
 plot(abs(V(1:length(t)/4)),'linewidth', 1)
 % approx -- no square root!
 approx_mag = max(abs(real(V)),abs(imag(V))) + 0.4* min(abs(real(V)),abs(imag(V)));
 hold on
 plot(approx_mag(1:length(t)/4),'r')
 legend('exact mag', 'approx mag')
 
 %compare exact and approx log
 % see:
 %  Generation of Products and Quotients Using 
 %  Approximate Binary Logarithms for Digital Filtering Applications
 %  IEEE Transactions on Computers 1970 vol.19 Issue No.02 - February
 logV = log2(abs(V));
 subplot(3,1,3);
 % plot the exact
 plot(logV(1:length(t)/4),'linewidth', 1)
 % get approx log of approx magnitude
 for i=1:length(t)
    sx = approx_mag(i); %abs(V(i));
    y=1; ly=0;
    while(sx>1)
        y=y*2; ly=ly+1; sx=floor(sx/2);
    end   
    % following approx +/- 0.04
    approx_log(i) = ly + (approx_mag(i)-y)./(y) + 0.043;
 end
 hold on
 %and plot it
 plot(approx_log(1:length(t)/4),'r')
 legend('exact log mag', 'approx log mag')
 
 
 