% 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