+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP Prime Software Library (/forum-15.html)
+--- Thread: Spectrum Analysis (Fourier) (/thread-4539.html)
Spectrum Analysis (Fourier) - toshk - 08-17-2015 04:51 PM
It works by plotting a function at F0.
Spect() to plot the spectra.
eg. SIN(X) at F0 plotted within a range of period T=2*pi.
PHP Code:
#pragma mode( separator(.,;) integer(h32) )
export Spect()
BEGIN
local dmin,dmax,dx,n;
local L0,L1,L4,L5;
L0:={}; L1:={}; L4:={}; L5:={};
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));
END;
L5:=makelist(X,X,0,(1/dx)/2,(1/dx)/1025);
L0:={}; L1:={};
startapp("Statistics 1Var");
FOR n FROM 1 TO size(L5) STEP 1 DO
L1[n]:=L4[n]; end;
D1:=L5; D2:=L1;
Statistics_1Var.H1:={'D2','D1','5'};
Statistics_1Var.CHECK(1);
startview(1,1);
END;
RE: Spectrum Analysis (Fourier) - toshk - 08-15-2018 07:35 PM
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;
[attachment=6218]
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;
#endRE: Spectrum Analysis (Fourier) - fabiofilippini - 03-18-2021 02:06 PM
Hi, how can I use your program with Ish?