Spectrum Analysis (Fourier)

[toshk]

compile compile compile compile!!!!
The plot you see after a successful run is Magnitude Plot;
For Phase plot; plot D3;

I
Spect()
Under this notation in F0 with Xrange of a period or integer multiple of T.
eg. F0 for sin(X) in Xrng 0 to 2*pi or -pi to +pi
eg. 1*Heaviside(X)+-1*Heaviside(X-1) Xrng 0 to 2; this defines PiecewiseFunction of X=1 between 0 and 1; X=0 elsewhere repeating at a T=2;

II
spect(lowerLimit, Period)
Under this notation plot in F0; don't care about Xrng but specify the range in function call.
eg. F0 for sin(X); call function spect(0,2*pi) OR spect(-pi,2*pi); bcos sinX T=2pi always.
the ~Fourier plot can be found under Function app. Press plot under Function App.
eg. 1*Heaviside(X)+-1*Heaviside(X-1); PiecewiseFunction of X=1 between 0 and 1; X=0 elsewhere;
function call spect(0,2) for T=2; check plot on Function App;
   

Code:


#pragma mode( separator(.,;) integer(h32) )
export Spect()
BEGIN
//Under this notation in F0 with Xrange of a period or integer multiple of T.
//eg F0 for sin(x) in Xrng 0 to 2*pi or -pi to +pi
//eg 1*Heaviside(X)+-1*Heaviside(X-1) Xrng 0 to 2; this defines PiecewiseFunction of X=1 between 0 and 1; X=0 elsewhere repeating at a T=2;
local dmin,dmax,dx,n;
local L0,L1,L4,L5,d3;
L0:={}; L1:={}; L4:={}; L5:={}; d3:={};
dmin:=Function.Xmin;
dmax:=Function.Xmax;
dx:=(dmax-dmin)/1024;
L0:=makelist(X,X,dmin,dmax,dx);
L1:=F0(L0);
L4:=mat2list(fft(L1));
FOR n FROM 1 TO size(L4) STEP 1 DO
L4[n]:=approx(ABS(L4[n])*(2/1024));
d3[SIZE(L4)]:=0;
IF TYPE(L4[n])==3 THEN d3[n]:=ARG(L4[n]); END;
END;
L5:=makelist(X,X,0,(1/dx)/2,(1/dx)/1025);
L0:={}; L1:={};
startapp("Statistics 1Var");
L1:=L4[{1,size(L5)}]; L1[1]:=L1[1]/2;
D1:=L5; D2:=L1;
Statistics_1Var.D3:=mat2list(d3[{1,size(L5)}]);
Statistics_1Var.H1:={'D2',"",'5'};
Statistics_1Var.CHECK(1);
startview(1,1);
END;

#cas
spect(ll,TT)
BEGIN
//spect(lowerLimit, Period)
//Under this notation plot in F0; don't care about Xrng but specify the range in function call.
//eg F0 for sin(X); call function spect(0,2*pi) OR spect(-pi,2*pi); bcos sinX T=2pi always.
//the ~Fourier plot can be found under Function app. Press plot under Function App.
//eg 1*Heaviside(X)+-1*Heaviside(X-1); PiecewiseFunction of X=1 between 0 and 1; X=0 elsewhere;
//function call spect(0,2) for T=2; check plot on Function App; 
 LOCAL aa,n,d3,d1,na;
 aa:={}; d3:={}; d1:={}; na:=0;

IF Q==4352 THEN aaa:=Function.F8(x); IF ll<0 THEN aaa:=Function.F8(x+ll); na:=ll; ll:=0; END;


ELSE aaa:=Function.F0(x); IF ll<0 THEN aaa:=Function.F0(x+ll); na:=ll; ll:=0; END; END;
  
 IFERR aa:=fourier_cn(aaa,'x',ABS(TT),{0,1,2,3,4,5,6,7,8,9},ll); THEN aa:={}; END;

 //IF SIZE(aa)==0 THEN L0:=MAKELIST(X,X,ll,TT,(TT-ll)/256); L1:=subst(Function.F0(x),x,L0); L2:=mat2list(fft(L1))/256; L2[1]:=L2[1]/2; END; 
 L2:=suppress(aa,1); aa:=prepend(mat2list(L2*2),approx(aa[1])); d3[SIZE(aa)]:=0; d1[SIZE(aa)]:=0; d1[1]:=aa[1];
 FOR n FROM 2 TO size(aa) STEP 1 DO
 IF RE(aa[n])<>0 AND IM(aa[n])<>0  THEN d3[n]:=IM(aa[n])/RE(aa[n]); d1[n]:=ABS(aa[n]); END;
 IF RE(aa[n])<>0 AND IM(aa[n])==0  THEN d3[n]:=pi/2; d1[n]:=RE(aa[n]); END;
 IF RE(aa[n])==0 AND IM(aa[n])<>0  THEN d3[n]:=0; d1[n]:=-IM(aa[n]); END;
 END;
 L2:=d1;
 //f8:=aa[1]+sum(d1[t]*SIN((2*π*(t-1)*'x'/TT)+d3[t]),t,2,SIZE(aa)); //write function in sinX only
 f9:=aa[1]+sum(-IM(aa[t])*SIN(2*π*(t-1)*'x'/TT)+RE(aa[t])*SIN((2*π*(t-1)*'x'/TT)+pi/2),t,2,SIZE(aa)); //write function in sinX+cosX
 //f10:=laplace(aa[1])+sum(-IM(aa[t])*laplace(SIN(2*π*(t-1)*'x'/TT))+RE(aa[t])*laplace(SIN((2*π*(t-1)*'x'/TT)+pi/2)),t,2,SIZE(aa));
 f10:=(aa[1]/'x')+sum(-IM(aa[t])/(2*π*(t-1)/TT)*(1)/('x'^2+1)+RE(aa[t])/((2*π*(t-1)/TT))*('x')/('x'^2+1),t,2,SIZE(aa));
 //f10:=(aa[1]/'x')+sum(-IM(aa[t])*(2*π*(t-1)/TT)/('x'^2+(2*π*(t-1)/TT)^2)+RE(aa[t])*('x')/('x'^2+(2*π*(t-1)/TT)^2),t,2,SIZE(aa));
 f9:=subst(f9,x,x-na);
startapp("Function"); 
 F9(x):=""; F9(x):=f9; 
 //F8(x):=""; F8(x):=f8;
 Function.CHECK(9);
 startapp("Statistics 1Var");
 Statistics_1Var.D2:=ABS(aa); d3:={0,0,0,0,0,0,0,0,0,0};
 FOR n FROM 2 TO size(aa) STEP 1 DO
 IF TYPE(aa[n])==3 AND IM(aa[n])>epsilon THEN d3[n]:=ARG(aa[n]); END;
 END;
 Statistics_1Var.D1:=mat2list(L2);
 Statistics_1Var.D3:=mat2list(d3);
 Statistics_1Var.H1:={"D2","",5};  Statistics_1Var.H2:={"D3","",5};
 Statistics_1Var.CHECK(1);
 startview(1,1);
 startview(10,1);
END;

//F8*F0 Convolution
Convd(ll1,TT1,ll2,TT2)
BEGIN
 LOCAL f1z,d33,d333;
 Q:=4352; 
 spect(ll1,TT1); d33:=f10; Q:=0; 
 spect(ll2,TT2); d333:=f10; Q:=0; 
 f9:=ilaplace(d33*d333*(1-e^(-x*TT1))*(1-e^(-x*TT2)),x,x);
 IF instring(string(f9),"ilaplace") THEN f1z:=Residues(d33*d333); f9:=ilaplace(f1z);END;
 startapp("Function");
 Xmin:=evalf(MIN(ll1,ll2)); Xmax:=2*evalf(gcd(exact(TT1),exact(TT2)));
 F9(x):=""; F9(x):=f9; Function.CHECK(9);
END;
#end