#! /bin/octave
% GNU Octave / MATLAB
% Benutzer:mik81 @ de.wikipedia.org



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



%
% Ermittle Signalform
%
 
function [square, realDutycycle] = waveform(dutycycle)
 
  global periods;
  global stepsPerPeriod;
    
  
  realDutycycle = 0;
  
  
  
  steps = periods*stepsPerPeriod;


  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
 
 
  dutycycleString = sprintf("%0.3f", dutycycle);
 
  plot((0:1:steps/(0.25*periods)-1)/stepsPerPeriod
      , square(1,1:steps/(0.25*periods)), [";Tastgrad: " , dutycycleString, ";"], "linewidth", 2);
  axis ([0, 0.25*periods, -0.2, +1.2]);
 
  axis ("nolabel");
  axis("labely", "ticxy");
%  xlabel ("");
  set(gca(), "linewidth", 1.5);


endfunction

%
% Berechne Spektrum
%


function [dft, envelop] = spektrum(square, realDutycycle)

  global periods;
  global stepsPerPeriod;
  global timebase;
  steps = periods*stepsPerPeriod;

  freqSteps = steps/2;
 
  %
  % Frequenz = freq*Samplerate*steps
  % n * Fenster
 
  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; 
 
  endfor
 
  dft = sqrt(dft);
 
  dft = dft/(periods*stepsPerPeriod);
  
  temp = 1/30; % gib nen sinnvollen Namen
 
  envelop = dft(1,1) * abs(sinc((0:temp:freqSteps)*realDutycycle));
  if(realDutycycle > 0.5)
    envelop2 = (1-dft(1,1)) * abs(sinc((0:temp:freqSteps)*(1-realDutycycle)));
  endif
 
  dftSelected = [dft(1,1)];
  for selected=periods+1:periods:freqSteps+1
    dftSelected = [dftSelected, dft(1,selected)];
  endfor
  
  dft = dftSelected;

 
  plot([0, (periods+1:periods:freqSteps+1)/periods-1/periods], dft(1, :), "or;Spektallinien;", "linewidth", 1.5);     
  hold("on");

  h = stem([0, (periods+1:periods:freqSteps+1)/periods-1/periods], dft(1, :), "r");
  set(h, "linewidth", 1.5);
  
  if(exist("envelop2", "var") != 1)
      plot((0:temp:freqSteps), envelop, ";Einhuellende;", "linewidth", 1.5);     

  else
    plot((0:temp:freqSteps), envelop, ";Einhuellende;", "linewidth", 0.5);     
    plot((0:temp:freqSteps), envelop2, "g;Invertierte Einhuellende;", "linewidth", 1.5);     
  endif

  hold("off");

  axis ([0, 20, 0, 1.1]);
  axis("labelxy", "ticxy", "on");
  xlabel ("Harmonische der Grundfrequenz");
  set(gca(), "linewidth", 1.5);
  %axis ("ticx");
  
endfunction

%
%
%  M  A  I  N
%
%

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

begin_cputime = cputime();

frames = 10; % Teilung
frames = 256; % 2^n

frames = frames-1; % Zaunpfahl

Dimensions = "-S1600,900";

global stepsPerPeriod = 512; % 2^n wegen FFT
global periods = 8; % 2^4
steps = periods*stepsPerPeriod;
global timebase = 0:1:steps-1;
size(timebase)
% dutycycle = 0.05;

mkdir("./png");
axis("labelxy", "ticxy", "on");
 
for count=1:1:frames % Hack: vermeide Tastgrad 1.00

  close();

  countString = sprintf("%03i", count);
  dutycycle = count/(frames+1);
 
 
  printf("\nCalculating waveform: %i\n", count);

  subplot(2,3,1, "align");
  [square, realDutycycle] = waveform(dutycycle);
   
%  print(["./png/waveform", countString, ".png"], "-dpng");

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

  subplot(1,3,2:3, "align");
  [dft, envelop] = spektrum(square, realDutycycle);

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


printf("Total cpu time: %f seconds\n", cputime()-begin_cputime);
printf("Finished\n");


% EOF