% GNU Octave / MATLAB
% Benutzer:mik81 @ de.wikipedia.org

%
%
% Große Frage:
%   Welche "amplitude" haben die Spektren?
%
%

printf("mik81 @ de.wikipedia.org\n\n");
printf("Spectrum of rectengular waveform\n\n");

stepsPerPeriod = 128; % 2^8
periods = 16; % 2^4
steps = periods*stepsPerPeriod;
timebase = 0:1:steps-1;
size(timebase)
% dutycycle = 0.05;

mkdir("./data");
axis("labelxy", "ticxy", "on");

for count=1:1:40

  countString = sprintf("%03i", count);
  dutycycle = count/80;
  dutycycleString = sprintf("%0.3f", dutycycle);
  realDutycycle = 0;

  %
  % set up waveform
  %

  printf("Calculating waveform: %i\n", count);

  for period=0:1:periods-1
    for step=0:1:stepsPerPeriod-1
      if ( step >= dutycycle *stepsPerPeriod)
        square = [square, 0.];
        if ( realDutycycle == 0)
          realDutycycle = step/stepsPerPeriod
          dutycycle
        endif
      else
        if (period == 0 && step == 0)
          square = [+1.];
        else
          square = [square, +1.];
        endif
      endif
    endfor
  endfor

  for step = 0:1:steps-1
    if ( step == 0)
      sinus = [sin(2*pi*0)];
    else
      sinus = [sinus, sin(2*pi*step/stepsPerPeriod)];
    endif
  endfor

  sinus = 0:1:steps-1;
  sinus = sin(2*pi*sinus/stepsPerPeriod);

  % size(sinus)
  % size(square)

  plot((0:1:steps/(0.25*periods)-1)/stepsPerPeriod
      , square(1,1:steps/(0.25*periods)), [";Tastgrad: " , dutycycleString, ";"]);
  axis ([0, 0.25*periods, -0.2, +1.2]);

  axis ("nolabel");
  axis("labely", "ticxy");
  xlabel ("");
  __gnuplot_set__ size 0.42,0.42;
  print(["./data/waveform", countString, ".png"], "-dpng");
  __gnuplot_set__ size 1.,1.;

  

  %
  % DFT
  %

  printf("Calculating DFT: %i\n", count);

  freqSteps = steps/2;

  %
  % Frequenz = freq*Samplerate*steps
  % n * Fenster

  % remove DC-part
  %square = square .- sum(square(1,:));

  for freq=0:1:freqSteps
    pad =  square .* sin(freq*2*pi*timebase/length(timebase));
    dft(1, freq+1) = sum( pad(1,:) );
    pad =  square .* cos(freq*2*pi*timebase/length(timebase));
    % cos^2 + sin^2
    dft( 1, freq+1) = dft(1, freq+1)^2 + sum( pad(1,:) )^2; 

  %  normalize = max(dft(1, freq+1), normalize);

    pad =  sinus .* sin(freq*2*pi*timebase/length(timebase));
    dftSignal(1, freq+1) = sum( pad(1,:) );
    pad =  sinus .* cos(freq*2*pi*timebase/length(timebase));
    % cos^2 + sin^2
    dftSignal(1, freq+1) = dftSignal( 1, freq+1)^2 + sum( pad(1,:) )^2;

  endfor

  dft = sqrt(dft);

  dft = dft/(periods*stepsPerPeriod);
  dftSignal = dftSignal/(periods*stepsPerPeriod);

  %envelop = dft(1,1) * abs(sinc((0:1/10:freqSteps)*dutycycle));
  envelop = dft(1,1) * abs(sinc((0:1/10:freqSteps)*realDutycycle));
  
%size(dft)

  %semilogx((0:1:steps-1)/stepsPerPeriod, square);
  %semilogx((1:1:freqSteps-1)/periods, dft(2:freqSteps,1));

  % plot

  plot((0:1:freqSteps)/periods, dft(1, 1:freqSteps+1), ["+1;Spektallinien;"]
      , (0:1:freqSteps)/periods, dft(1, 1:freqSteps+1), ["^1;;"]
      , (0:1/10:freqSteps), envelop, ";Einhuellende;");

  dftSelected = [dft(1,1)];
  for selected=periods+1:periods:freqSteps+1
    dftSelected = [dftSelected, dft(1,selected)];
  endfor

  %size(dftSelected)
  %size([0, (periods+1:periods:freqSteps+1)/periods])
  plot([0, (periods+1:periods:freqSteps+1)/periods-1/periods], dftSelected(1, :), ["+1;Spektallinien;"]
      ,[0, (periods+1:periods:freqSteps+1)/periods-1/periods], dftSelected(1, :), ["^1;;"]
      , (0:1/10:freqSteps), envelop, ";Einhuellende;");

  %semilogx((1:1:freqSteps-1)/periods, dftSignal(2:freqSteps,1));
  %plot((1:1:freqSteps/8-1)/periods, dftSignal(2:freqSteps/8,1));
  axis ([0, 20, 0, 0.6]);
  axis("labelxy", "ticxy", "on");
  xlabel ("Vielfache der Grundfrequenz");
  %axis ("ticx");

 
  print(["./data/spectrum", countString, ".png"], "-dpng");

  %if ( freq == 8)
  %    fetch = pad;
  %    sum(pad(1,:))
  %    dft(1, freq+1)
  %  endif
  %plot(0:1:length(fetch)-1, fetch, [";fetch", dutycycleString ,";"] );
  %print(["./data/fetch", countString, ".png"], "-dpng");

endfor