function geortzel_test

  f_sample = 8192;
  f_target = 2048;
%  f_target = 2000.0;

  N = 200;
  
  largest_mag = 0;

%t = 1/f_sample:(1/f_sample):N/f_sample;

%x = sin(2*pi*t*f_signal);
%close all;
%plot(t,x);

  h = hamming_coeff(N);
  %%%plot(h);pause
  k = N * f_target / f_sample;
  realW = cos(2.0*pi*k/N);
  imagW = sin(2.0*pi*k/N);
  
  fixed_one = 65536;
  log2_one_rep = int32(log2((fixed_one)));
  
  h8 = int32(h*fixed_one);
  
  %[fid,msg] = fopen('C:\files\Projects\msp430-workspace\Goertzel\coefficients.h','w');
  [fid,msg] = fopen('C:\files\Projects\msp430-goerzel\coefficients.h','w');
  fid
  msg
  fprintf(fid,'// center frequency: %dHz\r\n',f_target);
  fprintf(fid,'// sample frequency: %dHz\r\n',f_sample);
  fprintf(fid,'const uint8  window_size = %d;\r\n',N);
  fprintf(fid,'const int32 realW = %d;\r\n',int32(realW*fixed_one));
  fprintf(fid,'const int32 imagW = %d;\r\n',int32(imagW*fixed_one));
  fprintf(fid,'const uint8  log2_one_rep = %d;\r\n',log2_one_rep);
  fprintf(fid,'const uint16 scaling[] = {\r\n');
  for n=1:length(h8)
    if (h8(n) == fixed_one)
      fprintf(fid,'  %d,\r\n',h8(n)-1);
    else
      fprintf(fid,'  %d,\r\n',h8(n));
    end
  end
  fprintf(fid,'};\r\n');
  fclose(fid);
  
  %plot(h8);
  %pause
  frequencies = 0:8:4000;
  for i=1:length(frequencies)
    f = frequencies(i);
    x = generate_signal(f_sample, f, N);
    
    x8 = int32(512 + ((x+1)/2) * 128);
    %plot(x8);
    %pause
    
    %%%h=ones(length(x));
    [zr,zi] = goertzel(x,h,realW,imagW);
    %amplitude(i) = abs(zr+sqrt(-1)*zi);

    % [zr,zi] = goertzel88(x8,h8,int32(realW*fixed_one),int32(imagW*fixed_one));
    amplitude(i) = abs(zr)+abs(zi);
    % amplitude(i) = goertzel_ratio(x,h,realW,imagW);
    %%%amplitude(i) = goertzel_ratio_fixed(x8,h8,int32(realW*fixed_one),int32(imagW*fixed_one));
    % amplitude(f) = zr^2 + zi^2;
    update_mag(amplitude(i));
    % [ f amplitude(i) ]
    % pause;
  end
  
  plot(frequencies,amplitude,'')
  pause
  largest_mag

  function x = generate_signal(f_sample, f_signal, N)
    t = 1/f_sample:(1/f_sample):N/f_sample;

    x = sin(2*pi*t*f_signal);
  end

  function r = goertzel_ratio(x,h,realW,imagW)
    N = length(x);
    
    d1 = 0.0;
    d2 = 0.0;
    a  = 0;
    for n=1:N
      s = h(n) * x(n);
      a = a + abs(s);
      y  = h(n) * x(n) + 2*realW*d1 - d2;
      update_mag(y);
      % Y(n) = y; % just for plotting.
      d2 = d1;
      d1 = y;
    end
    
    resultr = realW*d1 - d2;
    resulti = imagW*d1;

    update_mag(resultr);
    update_mag(resulti);
        
    r = (abs(resultr)+abs(resulti))/a
    
    % plot(Y);
  end

  function [resultr,resulti] = goertzel(x,h,realW,imagW)
    N = length(x);
    
    d1 = 0.0;
    d2 = 0.0;
    for n=1:N
      y  = h(n) * x(n) + 2*realW*d1 - d2;
      update_mag(y);
      % Y(n) = y; % just for plotting.
      d2 = d1;
      d1 = y;
    end
    
    resultr = realW*d1 - d2;
    resulti = imagW*d1;
    
    update_mag(resultr);
    update_mag(resulti);
    % plot(Y);
  end

  function r = goertzel_ratio_fixed(x,h,realW,imagW)
    N = length(x);
    
    d1 = int32(0);
    d2 = int32(0);
    a  = int32(0);
    for n=1:N
      %[n h(n) x(n) h(n)*x(n) (h(n)*x(n)/fixed_one) ]
      %pause;
      y1 = h(n) * (x(n) - 512);
      update_mag(y1);
      y1 = y1 / fixed_one;

      a = a + abs(y1);
      
      y2 = realW * d1;
      update_mag(y2);
      y2 = y2 / (fixed_one/2);

      y = y1 + y2 - d2;
      update_mag(y);
      % y  = h(n) * x(n) + 2*realW*d1 - d2;
      d2 = d1;
      d1 = y;
    end
    
    resultr = (int32(realW)*d1 / (fixed_one)) - d2;
    resulti = int32(imagW)*d1 / (fixed_one);
    update_mag(resultr);
    update_mag(resulti);
    
    [abs(resultr)+abs(resulti) a]
    
    r = (abs(resultr)+abs(resulti))*fixed_one/a;
    r = double(r) / fixed_one;

    %plot(Y);
    %pause
  end

  % 8.8 fixed-point
  function [resultr,resulti] = goertzel88(x,h,realW,imagW)
    N = length(x);
    
    d1 = int32(0);
    d2 = int32(0);
    for n=1:N
      %[n h(n) x(n) h(n)*x(n) (h(n)*x(n)/fixed_one) ]
      %pause;
      y1 = h(n) * x(n);
      update_mag(y1);
      y1 = y1 / fixed_one;

      y2 = realW * d1;
      update_mag(y2);
      y2 = y2 / (fixed_one/2);

      y = y1 + y2 - d2;
      update_mag(y);
      % y  = h(n) * x(n) + 2*realW*d1 - d2;
      Y(n) = y; % just for plotting.
      d2 = d1;
      d1 = y;
    end
    
    resultr = (int32(realW)*d1 / fixed_one) - d2;
    resulti = int32(imagW)*d1 / fixed_one;
    update_mag(resultr);
    update_mag(resulti);

    %plot(Y);
    %pause
  end


  function [resultr,resulti] = geortzel8(x8,h8,f_sample,f_target)
    N = length(x);
    k = N * f_target / f_sample;
    half_realW = cos(2.0*pi*k/N);
    imagW = sin(2.0*pi*k/N);
    
    d1 = 0.0;
    d2 = 0.0;
    for n=1:N
      y  = h(n) * x(n) + 2*half_realW*d1 - d2;
      update_mag(y);
      % Y(n) = y; % just for plotting.
      d2 = d1;
      d1 = y;
    end
    
    resultr = half_realW*d1 - d2
    resulti = imagW*d1
    
    update_mag(resultr);
    update_mag(resulti);
    % plot(Y);
  end

  function h = hamming_coeff(N)
    for n=0:N-1
      h(n+1) = 0.54 - 0.46 * cos(2*pi*n/N);
    end
  end

  function update_mag(x)
    if (x == 32767)
      x
    end
    largest_mag = max(largest_mag,x);
  end

end