Classic Fourier Series

[salvomic]

hi,
I need something to calculate, in HP Prime, the n-th Fourier coefficient (like in HP50: fourier(x^2,0) -> 4π^2 / 3).
I can't find it in catalog of Prime. There are only FFT and its inverse...

help :-)

Salvo

[Snorre]

Hello,

why not define it by yourself:
a(f,k):=int(f(t)*cos(k*t),t=-π..π)/π
b(f,k):=int(f(t)*sin(k*t),t=-π..π)/π
or something similiar.
(Beware: this may suffer from a little issue which will be fixed in next software version.)

Greetings

[salvomic]

(01-30-2015 05:20 PM)Snorre Wrote:  Hello,

why not define it by yourself:
a(f,k):=int(f(t)*cos(k*t),t=-π..π)/π
b(f,k):=int(f(t)*sin(k*t),t=-π..π)/π
or something similiar.
(Beware: this may suffer from a little issue which will be fixed in next software version.)

Greetings

thank you, I'll think about this suggestion...
I hope also the little issue could solved soon

However, it would be likely to have such a function included in the CAS, if it's possible ;-)

cheers
Salvo

[salvomic]

define the functions it works, but I get first a lot of warnings, like those claiming the use of Taylor or subs() and so on...
ex.: "Warning, evalueting univariate expression(value) like if expression was a function. You should write subst(t^2,t,taylorx21" if I try with t^2 as function

[Han]

(01-30-2015 06:48 PM)salvomic Wrote:  define the functions it works, but I get first a lot of warnings, like those claiming the use of Taylor or subs() and so on...
ex.: "Warning, evalueting univariate expression(value) like if expression was a function. You should write subst(t^2,t,taylorx21" if I try with t^2 as function

This may have to do with using an expression as a function. Are you passing a CAS function as your argument? Or are you passing a CAS expression?

[salvomic]

(01-30-2015 07:23 PM)Han Wrote:  This may have to do with using an expression as a function. Are you passing a CAS function as your argument? Or are you passing a CAS expression?

I tried this a(t^2,1)

***
then some variables were a mess, until I purge some of them: d, g, v, used in another programs...

[Han]

(01-30-2015 07:32 PM)salvomic Wrote:  

(01-30-2015 07:23 PM)Han Wrote:  This may have to do with using an expression as a function. Are you passing a CAS function as your argument? Or are you passing a CAS expression?

I tried this a(t^2,1)

***
then some variables were a mess, until I purge some of them: d, g, v, used in another programs...

If you look at how a(f,k) was defined, the f is assumed to be a CAS function (notice the f(t) reference in how a(f,k) is defined). This means you need to pass a CAS function.

f(t):=t^2;
a(f,1);

The CAS allows one to use an expression as a function. (This does not mean that it is good practice.) That is, if we define f as an expression f:=t^2, and then do f(2), it will not only give a warning about using functional notation to evaluate an expression, but also give you the result 4.

[salvomic]

(01-30-2015 07:54 PM)Han Wrote:  If you look at how a(f,k) was defined, the f is assumed to be a CAS function (notice the f(t) reference in how a(f,k) is defined). This means you need to pass a CAS function.

f(t):=t^2;
a(f,1);

in this way it works (I used also variant with Int(...,0, 2π) ), but I tried this:
create a0 (first term), ak, bk variable with expressions
f(t):=|sin(t)|
I got warning about the absolute value...
then I got value for a0 (as a0:=(1/2π)Int(f(t),t,0,2π)
After this value for a0, ak, bk and also f=f(t) are reset and I must rewrite all three functions...
There is a reason? Is is possible to maintain the function

Perhaps is better a little program that give a list with (a0, ak, bk)?

[Han]

(01-30-2015 09:06 PM)salvomic Wrote:  

(01-30-2015 07:54 PM)Han Wrote:  If you look at how a(f,k) was defined, the f is assumed to be a CAS function (notice the f(t) reference in how a(f,k) is defined). This means you need to pass a CAS function.

f(t):=t^2;
a(f,1);

in this way it works (I used also variant with Int(...,0, 2π) ), but I tried this:
create a0 (first term), ak, bk variable with expressions
f(t):=|sin(t)|
I got warning about the absolute value...
then I got value for a0 (as a0:=(1/2π)Int(f(t),t,0,2π)
After this value for a0, ak, bk and also f=f(t) are reset and I must rewrite all three functions...
There is a reason? Is is possible to maintain the function

Perhaps is better a little program that give a list with (a0, ak, bk)?

I'm not sure I understand what you did; perhaps you could copy/paste the actual commands you used? I am not having any trouble with respect to things "disappearing"

[salvomic]

(01-31-2015 12:05 AM)Han Wrote:  I'm not sure I understand what you did; perhaps you could copy/paste the actual commands you used? I am not having any trouble with respect to things "disappearing"

ok, I'll try to copy/paste as soon as possible later, however I tried to do like you did. The issue was after the message "warning"...
Later I'll try again.

For now I'm trying also with a program like this:

Code:


#cas
fourcoeff(fun,k)
BEGIN
local a0:=(int(fun,t,0,2*PI))/(2*PI);
local ak:=(int(fun*cos(k*t),t,0,2*PI))/PI;
local bk:=(int(fun*sin(k*t),t,0,2*PI))/PI;
return a0,ak,bk;
END;
#end

it should give a0 (the first), ak and bk coefficient for k and a function.
It could be made better; if I miss the local command for variables it give formulas with something like 256 or 1024 added, I don't understand why...

Thank you for your effort!

Salvo

[Snorre]

Hello Salvo,

a program is always a good idea ;-)

Just two tips:

1st:
You should take care of the argument types your program expects.
Your code seems to expect fun to be an expression (not a function) with t (and only t) being the independent variable, so the typical usage looks like:
fourcoeff(t²,1) or
f(x):=x² and fourcoeff(f(t),1).
By that you have to always remember that fun must be an expression of t. And furthermore, t must not be defined already.
If you change your program slightly by

Code:

...
a0:=int(fun(t),t,0,2*PI)/(2*PI);
ak:=int(fun(t)*cos(k*t),t,0,2*PI)/PI;
bk:=int(fun(t)*sin(k*t),t,0,2*PI)/PI;
...

it'll expect fun to be a function (not an expression) of arity 1. The advantage of this approach is: it doesn't matter how you name the independent variable since your program will replace it by the local t anyway which is totally independent of any global t already defined. That is, you could do something like: t:=["foo","bar"] and g(x):=x² and fourcoeff(g,1).

2nd:
Your program returns a sequence (an object similiar to those created by the builtin seq function), not a vector/list.
That is totally OK as long as it was your intention. Just be aware that sequences behave quite differently when used in further calculations, for example 0+fourcoeff(...) vs. 0+[fourcoeff(...)].
If you want your program to return a vector or list, put the last expression return a0,ak,bk; within brackets or braces: return [a0,ak,bk]; or return {a0,ak,bk};.

Greetings

[Tim Wessman]

(01-30-2015 06:07 PM)salvomic Wrote:  However, it would be likely to have such a function included in the CAS, if it's possible ;-)

XCAS has 3 commands for this. When we was going through the list of commands at some point in the past to remove things that we did not have time to fully test or made no sense on the handheld calculator I suspect this were missed unintentionally.

Those commands are fourier_an, fourier_bn, fourier_cn defined as:

Code:

Returns the n-th Fourier coefficient an=2/T*integrate(f(x)*cos(2*pi*n*x/T),a,a+T).

Returns the n-th Fourier coefficient bn=2/T*integrate(f(x)*sin(2*pi*n*x/T),a,a+T).

Returns the n-th Fourier coefficient cn=1/T*integrate(f(x)*exp(-2*i*pi*n*x/T),a,a+T).

I assume those will do what you are looking for?

To everyone else, *please* point out missing commands from xcas if you think there is something that should be in Prime that is missing. Note that each command that is added however introduces additional overhead in the form of testing, documentation and so forth. Simply "turning everything on" does not make sense because of time contraints, the large numbers of synonyms/intermingled french command names, things that do not have compatible/valid syntax on the HP side, or things that just can't be supported for one reason or another at the moment (like writing out a wav file... yes, the cas will do that with the "writewav" command...)

[salvomic]

(01-31-2015 09:29 PM)Tim Wessman Wrote:  XCAS has 3 commands for this. When we was going through the list of commands at some point in the past ...

Those commands are fourier_an, fourier_bn, fourier_cn ...

I assume those will do what you are looking for?

yes, Tim!
I mean those commands, all three (cn is for the exponential series)...
So, they will be no more in the CAS, never in the next update?
Ok, then I must just my little program (see above), also if it would be better to have them integrated...

For me there few things that I could see missed in confront of my "old" HP 50g: these ones and the command "image" that in the HP50g is in the section "Linear Application" (with mkisom, ker...), for a matrix...
Another good thing to have would be "Directional Derivative" (see my post here)
For the rest Prime has already almost everything (excluding Graph3D and QPI, for example)
I'm adding only something to calculate inverse (pseudo inverse) of a matrix, and few other useful programs...

thank you,
Salvo

[salvomic]

(01-31-2015 09:11 PM)Snorre Wrote:  Hello Salvo,

a program is always a good idea ;-)

Just two tips:

...

hi friend,
your tips are very welcome!
I'm changing my program in this way.
I'm still a newcomer in this programming language, so I love to learn more and more tips!

Yes, however, I would export a list of three arguments: a0 (however this one is always the same), ak, bk...

Last: doing so:

Code:


a0:=int(fun(t),t,0,2*PI)/(2*PI);
ak:=int(fun(t)*cos(k*t),t,0,2*PI)/PI;
bk:=int(fun(t)*sin(k*t),t,0,2*PI)/PI;

How can we avoid this error?
I set g:=x^2 then input fourcoeff(g(x),1) or fourcoeff(x^2,1) and I get "warning, evaluating invariable expression(value) like if expression was a function. You should write subst(x^2,x,t)"...
Sure I miss something else
I would better have the possibility to input both functions or expressions :-)

For the rest, I would also use the program (or the function we discussed above) with some "piecewise" functions.
So I tried to define a function for "square wave", f(t):=PIECEWISE(-π<t<0,-1,t=0,0,t=-π,0,t=π,0,0<t<π,1); this function seem to be ok, but if I use it with the program or with other integrate function I get "No case applies"...
Is this the way to use "Piecewise" function or I miss something?

Thank you!

Salvo

[Han]

(01-31-2015 09:53 PM)salvomic Wrote:  Yes, however, I would export a list of three arguments: a0 (however this one is always the same), ak, bk...

Last: doing so:

Code:


a0:=int(fun(t),t,0,2*PI)/(2*PI);
ak:=int(fun(t)*cos(k*t),t,0,2*PI)/PI;
bk:=int(fun(t)*sin(k*t),t,0,2*PI)/PI;

How can we avoid this error?
I set g:=x^2 then input fourcoeff(g(x),1) or fourcoeff(x^2,1) and I get "warning, evaluating invariable expression(value) like if expression was a function. You should write subst(x^2,x,t)"...
Sure I miss something else
I would better have the possibility to input both functions or expressions :-)

Again, you need to understand the difference between a function and an expression.

g:=x^2; <-- this creates an expression whose name is g
g(x):=x^2; <-- this creates a function whose name is g

a0:=int(fun(t),t,0,2*PI)/(2*PI); <-- this creates a variable named a0 whose value is the result of integrating a function named 'fun' (which hopefully you've already pre-defined); if 'fun(t)' is not defined, then a0 is an expression (otherwise a0 is presumably a numerical value)

ak:=int(fun(t)*cos(k*t),t,0,2*PI)/PI; <-- this creates a variable named ak whose value is the result of integrating the product of the function 'fun(t)' and 'cos(k*t)'; if 'fun' and 'k' were predefined prior to creating ak, then ak is just a numerical value; otherwise ak is an expression (same as a0)

Neither a0 nor ak are functions when defined this way. However, the CAS will let you get away with using expressions as if they were functions -- it will just always throw warnings/complaints about using expressions as if they were functions. So don't do it if you don't want to see the warnings.

ak(f,k):=int(f(t)*cos(k*t),t,0,2*PI)/(2*PI); <-- this creates a function named ak, whose arguments are a CAS function f and a constant k; that means 'f' must be a function (see the g(x):=x^2 example), and NOT an expression (see the g:=x^2 example)

Once you understand the difference between the two, your warnings will go away.

[rprosperi]

(02-01-2015 03:10 AM)Han Wrote:  Again, you need to understand the difference between a function and an expression.

g:=x^2; <-- this creates an expression whose name is g
g(x):=x^2; <-- this creates a function whose name is g

a0:=int(fun(t),t,0,2*PI)/(2*PI); <-- this creates a variable named a0 whose value is the result of integrating a function named 'fun' (which hopefully you've already pre-defined); if 'fun(t)' is not defined, then a0 is an expression (otherwise a0 is presumably a numerical value)

ak:=int(fun(t)*cos(k*t),t,0,2*PI)/PI; <-- this creates a variable named ak whose value is the result of integrating the product of the function 'fun(t)' and 'cos(k*t)'; if 'fun' and 'k' were predefined prior to creating ak, then ak is just a numerical value; otherwise ak is an expression (same as a0)

Neither a0 nor ak are functions when defined this way. However, the CAS will let you get away with using expressions as if they were functions -- it will just always throw warnings/complaints about using expressions as if they were functions. So don't do it if you don't want to see the warnings.

ak(f,k):=int(f(t)*cos(k*t),t,0,2*PI)/(2*PI); <-- this creates a function named a0, whose arguments are a CAS function f and a constant k; that means 'f' must be a function (see the g(x):=x^2 example), and NOT an expression (see the g:=x^2 example)

Once you understand the difference between the two, your warnings will go away.

Han - thanks for these explanations/clarifications; your detailed examples are very useful for learing to use CAS, and some of the more exotic features.

fyi - For your final example, in the explanation 'a0' should be 'ak', that is if I understood it correctly.

[Han]

(02-01-2015 03:22 AM)rprosperi Wrote:  fyi - For your final example, in the explanation 'a0' should be 'ak', that is if I understood it correctly.

You are quite right; I've gone back and made the proper edits. Thank you!

[salvomic]

(02-01-2015 03:10 AM)Han Wrote:  ...
ak(f,k):=int(f(t)*cos(k*t),t,0,2*PI)/(2*PI); <-- this creates a function named ak, whose arguments are a CAS function f and a constant k; that means 'f' must be a function (see the g(x):=x^2 example), and NOT an expression (see the g:=x^2 example)

Once you understand the difference between the two, your warnings will go away.

thank you, you are very kind about this explanation, very clear, now, for me...
I'll try some variations in the program.

Have a nice Sunday!

[salvomic]

(01-31-2015 09:11 PM)Snorre Wrote:  a program is always a good idea ;-)

Just two tips:

Hi Snorre,
after some tries now I'm working on this idea, just for fun (and to understand):

Code:


#cas
fourcoeff(fun,k)
BEGIN
local nstr := instring(string(fun),"t");
// check if t is in string fun
local tipo:=type(fun);
// check the type: DOM_FUNC, DOM_SYMBOLIC, ...
local a0, ak, bk;
CASE
IF tipo = DOM_FUNC THEN
a0:=(int(fun(t),t,0,2*PI))/(2*PI);
ak:=(int(fun(t)*cos(k*t),t,0,2*PI))/PI;
bk:=(int(fun(t)*sin(k*t),t,0,2*PI))/PI;
return {a0,ak,bk};  // a list
END;
IF nstr != 0 THEN
// not t
a0:=(int(fun,t,0,2*PI))/(2*PI);
ak:=(int(fun*cos(k*t),t,0,2*PI))/PI;
bk:=(int(fun*sin(k*t),t,0,2*PI))/PI;
return {a0,ak,bk};  // a list
END;
DEFAULT return "Input must be function or expr with t";
END;
END;
#end

With this, after having defined g(t):=t^2, I get:
fourcoeff(g,1) -> ok, correct value
fourcoeff(t^2,1) -> ok, correct value, without any warning
fourcoeff(x^2,1) -> exit with the message.

This is a workaround, sure, and the the program could be strongly improved!

So the "type" control is working (it now handle a function in t and and expression in t), but it would better to set another control to say "if tipo=DOM_SYMBOLIC and the letter is not "t" substitute any letter like x, y, z... into "t") or reformat fun in fun(t) in some way but a efficient one (now the expression could have "t" without to be correct)...

Thank you for patience

[Snorre]

Hi,

another approach:

Code:

#cas
fourier(args):=
BEGIN
  LOCAL argv,argc,f,k;
  argv:=[args];
  argc:=SIZE(argv);
  f:=argv(1);
  k:=argv(argc);
  IF argc=3 THEN f:=zip('unapply',f,argv(2)); END;
  RETURN int(f(t)*e^(-i*k*t),t,-π,π)/π;  // <-- replace that by your a0/ak/bk formulas
END;
#end

This doesn't check types, but number of args, so usage is either fourier(expr,var,k) or fourier(func,k).