% test_fft2.m  test various Fourier series
function test_fft2
  clear
  format compact
  diary 'test_fft2_m.out'
  
  Pi=3.1415926535897932384626433832795028841971;
  n=16;
  rdat=zeros(n);
  idat=zeros(n);

  disp('test_fft2.m')
  for i=1:n % sin 2 th  other coefficients zero
    t=i-1.0;
    th = (2.0*2.0*Pi)*t/n;
    rdat(i) = sin(th);
    idat(i) = 0.0;
  end
  disp('sin 2 th data')
  cxdat=rdat+j*idat;
  for i=1:n 
    disp([num2str(i-1) ' real = ' num2str(real(cxdat(i))) ...
                      ', imag = ' num2str(imag(cxdat(i)))])
  end
  rslt=fft(cxdat);
  disp('sin 2 th FFT')
  for i=1:n
    disp([num2str(i-1) ' cos = ' num2str(real(rslt(i))) ...
                      ', sin = ' num2str(imag(rslt(i)))])
  end
  disp('antialias and normalize sin 2 th FFT')
  for i=2:n/2
    rslt(i)=conj((rslt(i)+conj(rslt(n+2-i)))/n);
  end
  rslt(1)=conj(rslt(1)/n);
  for i=1:n/2
    disp([num2str(i-1) ' cos = ' num2str(real(rslt(i))) ...
                      ', sin = ' num2str(imag(rslt(i)))])
  end
  disp(' ')
  disp(' ')

  for i=1:n % cos 2 th 
    t=i-1.0;
    th = (2.0*2.0*Pi)*t/n;
    rdat(i) = cos(th);
    idat(i) = 0.0;
  end
  disp('cos 2 th data')
  cxdat=rdat+j*idat;
  for i=1:n 
    disp([num2str(i-1) ' real = ' num2str(real(cxdat(i))) ...
                      ', imag = ' num2str(imag(cxdat(i)))])
  end
  rslt=fft(cxdat);
  disp('cos 2 th FFT')
  for i=1:n
    disp([num2str(i-1) ' cos = ' num2str(real(rslt(i))) ...
                      ', sin = ' num2str(imag(rslt(i)))])
  end
  disp('antialias and normalize cos 2 th FFT')
  for i=2:n/2
    rslt(i)=conj((rslt(i)+conj(rslt(n+2-i)))/n);
  end
  rslt(1)=conj(rslt(1)/n);
  for i=1:n/2
    disp([num2str(i-1) ' cos = ' num2str(real(rslt(i))) ...
                      ', sin = ' num2str(imag(rslt(i)))])
  end
  disp(' ')
  disp(' ')

  
  for i=1:n % 1.0 + cos 2 th 
    t=i-1.0;
    th = (2.0*2.0*Pi)*t/n;
    rdat(i) = 1.0+cos(th);
    idat(i) = 0.0;
  end
  disp('1.0 + cos 2 th data')
  cxdat=rdat+j*idat;
  for i=1:n 
    disp([num2str(i-1) ' real = ' num2str(real(cxdat(i))) ...
                      ', imag = ' num2str(imag(cxdat(i)))])
  end
  rslt=fft(cxdat);
  disp('1.0 + cos 2 th FFT')
  for i=1:n
    disp([num2str(i-1) ' cos = ' num2str(real(rslt(i))) ...
                      ', sin = ' num2str(imag(rslt(i)))])
  end
  disp('antialias and normalize 1.0 + cos 2 th FFT')
  for i=2:n/2
    rslt(i)=conj((rslt(i)+conj(rslt(n+2-i)))/n);
  end
  rslt(1)=conj(rslt(1)/n);
  for i=1:n/2
    disp([num2str(i-1) ' cos = ' num2str(real(rslt(i))) ...
                      ', sin = ' num2str(imag(rslt(i)))])
  end
  disp(' ')
  disp(' ')


  for i=1:n % square wave 4/Pi(cos th - 1/3 cos 3 th + 
            %  1/5 cos 5 th - 1/7 cos 7 th 
    t=i-1.0;
    v = 1.0;
    if t>=n/4 && t<3*n/4
      v = -1.0;
    end
    if t==n/4 || t==3*n/4
      v = 0.0;
    end
    rdat(i) = v;
    idat(i) = 0.0;
  end
  disp('square wave data, cos series')
  cxdat=rdat+j*idat;
  for i=1:n 
    disp([num2str(i-1) ' real = ' num2str(real(cxdat(i))) ...
                      ', imag = ' num2str(imag(cxdat(i)))])
  end

  disp('square wave generate, using cos')
  for i=1:n
    t=i-1.0;
    th = 2.0*Pi*t/n;
    v = (4.0/Pi)*(cos(th) - cos(3.0*th)/3.0 + cos(5.0*th)/5.0 - ...
        cos(7.0*th)/7.0 + cos(9.0*th)/9.0 - cos(11.0*th)/11.0 + ...
        cos(13.0*th)/13.0 - cos(15.0*th)/15.0 + cos(17.0*th)/17.0);
    disp([num2str(t) ' real = ' num2str(v)])
  end
  cxdat=rdat+j*idat;
  for i=1:n 
    disp([num2str(i-1) ' real = ' num2str(real(cxdat(i))) ...
                      ', imag = ' num2str(imag(cxdat(i)))])
  end
  rslt=fft(cxdat);
  disp('square wave FFT cos series')
  for i=1:n
    disp([num2str(i-1) ' cos = ' num2str(real(rslt(i))) ...
                      ', sin = ' num2str(imag(rslt(i)))])
  end
  disp('antialias and normalize square wave FFT cos')
  arslt=rslt;
  for i=2:n/2
    arslt(i)=conj((arslt(i)+conj(arslt(n+2-i)))/n);
  end
  arslt(1)=conj(arslt(1)/n);
  for i=1:n/2
    disp([num2str(i-1) ' cos = ' num2str(real(arslt(i))) ...
                      ', sin = ' num2str(imag(arslt(i)))])
  end
  disp(' ')
  disp(' ')

  cxdat=ifft(rslt);
  disp('square wave inverse FFT')
  for i=1:n 
    disp([num2str(i-1) ' real = ' num2str(real(cxdat(i))) ...
                      ', imag = ' num2str(imag(cxdat(i)))])
  end
  disp(' ')
  disp(' ')
  

  for i=1:n % square wave 4/Pi(sin th - 1/3 sin 3 th + 
            %  1/5 sin 5 th - 1/7 sin 7 th 
    t=i-1.0;
    v = 1.0;
    if t>n/2
      v = -1.0;
    end
    if t==n/2 || t==0
      v = 0.0;
    end
    rdat(i) = v;
    idat(i) = 0.0;
  end
  disp('square wave data, sin series')
  cxdat=rdat+j*idat;
  for i=1:n 
    disp([num2str(i-1) ' real = ' num2str(real(cxdat(i))) ...
                      ', imag = ' num2str(imag(cxdat(i)))])
  end

  disp('square wave generate, using sin')
  for i=1:n
    t=i-1.0;
    th = 2.0*Pi*t/n;
    v = (4.0/Pi)*(sin(th) + sin(3.0*th)/3.0 + sin(5.0*th)/5.0 + ...
        sin(7.0*th)/7.0 + sin(9.0*th)/9.0 + sin(11.0*th)/11.0 + ...
        sin(13.0*th)/13.0 + sin(15.0*th)/15.0 + sin(17.0*th)/17.0);
    disp([num2str(t) ' real = ' num2str(v)])
  end
  cxdat=rdat+j*idat;
  for i=1:n 
    disp([num2str(i-1) ' real = ' num2str(real(cxdat(i))) ...
                      ', imag = ' num2str(imag(cxdat(i)))])
  end
  rslt=fft(cxdat);
  disp('square wave FFT sin series')
  for i=1:n
    disp([num2str(i-1) ' cos = ' num2str(real(rslt(i))) ...
                      ', sin = ' num2str(imag(rslt(i)))])
  end
  disp('antialias and normalize square wave FFT sin')
  arslt=rslt;
  for i=2:n/2
    arslt(i)=conj((arslt(i)+conj(arslt(n+2-i)))/n);
  end
  arslt(1)=conj(arslt(1)/n);
  for i=1:n/2
    disp([num2str(i-1) ' cos = ' num2str(real(arslt(i))) ...
                      ', sin = ' num2str(imag(arslt(i)))])
  end
  disp(' ')
  disp(' ')
  

  for i=1:n % square wave 4/Pi(sin th - 1/3 sin 3 th + 
            %  1/5 sin 5 th - 1/7 sin 7 th 
    v = 1.0;
    if i>n/2
      v = -1.0;
    end
    rdat(i) = v;
    idat(i) = 0.0;
  end
  disp('square wave data, sin no zeros')
  cxdat=rdat+j*idat;
  for i=1:n 
    disp([num2str(i-1) ' real = ' num2str(real(cxdat(i))) ...
                      ', imag = ' num2str(imag(cxdat(i)))])
  end
  rslt=fft(cxdat);
  disp('square wave FFT sin no zeros')
  for i=1:n
    disp([num2str(i-1) ' cos = ' num2str(real(rslt(i))) ...
                      ', sin = ' num2str(imag(rslt(i)))])
  end
  disp('antialias and normalize square wave FFT sin')
  arslt=rslt;
  for i=2:n/2
    arslt(i)=conj((arslt(i)+conj(arslt(n+2-i)))/n);
  end
  arslt(1)=conj(arslt(1)/n);
  for i=1:n/2
    disp([num2str(i-1) ' cos = ' num2str(real(arslt(i))) ...
                      ', sin = ' num2str(imag(arslt(i)))])
  end
  disp(' ')
  disp(' ')
  diary off

end % test_fft2.m