#!/usr/bin/octave -qf

f1=1000;
f2=2000;
fmod=200;

N=2000;
t=(0:N)/N*4/fmod;

fsk=sin(2*pi*f1.*t).*(sin(fmod*2*pi.*t) > 0) \
    + sin(2*pi*f2.*t).*(sin(fmod*2*pi.*t)<0);

plot(t,fsk);
pause;

p1_sin=fsk.*sin(2*pi*f1.*t);
p1_cos=fsk.*cos(2*pi*f1.*t);
p2_sin=fsk.*sin(2*pi*f2.*t);
p2_cos=fsk.*cos(2*pi*f2.*t);

t1_win=t < 1/(f1); # Rechteckfenster
t2_win=t < 1/(f2);
## besser für ganzzahlige Frequenzverhältnisse:
#t2_win=t < 1/(f1);

out1=conv(p1_sin,t1_win).^2 + conv(p1_cos,t1_win).^2;
out2=(conv(p2_sin,t2_win).^2 + conv(p2_cos,t2_win).^2);
out2=out2.*(f2/f1)^2;  # Skalierung anpassen für t1_win <> t2_win

plot(t,out1(1:N+1),"",t,out2(1:N+1),"");
 
pause;