Search code examples
vectorizationcomplex-numbersscilabarray-broadcasting

Running feval() on a complex matrix


I am following Varodom Toochinda (a.k.a Dew)'s tutorial here. But I want to vectorize this piece of code:

f = linspace(0.01, 1, 1000); // frequency vector (Hz)
w = 2 * %pi * f;  // convert to rad/s
Pmag = zeros(1, 1000); Pph = zeros(1, 1000);
for k=1:1000, // compute P(jw) at each frequency point
  P = 1 / (-10 * w(k)^2 + 0.1 * %i * w(k)); 
  [Pmag(k), Pph(k)] = polar(P);
end

to

f = linspace(0.01, 1, 1000); // frequency vector (Hz)
w = 2 * %pi * f;  // convert to rad/s
P = 1 ./ (-10 * w.^2 + 0.1 * %i * w);
[Pmag, Pph_rad] = feval(P, polar);

but I get the

feval: Wrong number of output argument(s): 1 expected.

error. and if I try just the feval(P, polar) command, I get the

feval: Wrong type for input argument #1 : A real matrix expected.

error message. I would appreciate it if you could help me know what is the problem and how I can resolve it.


Solution

  • As (quite poorly but anyway) documented in the feval() help page :

    1. f is an [external] (function or routine) accepting on one or two arguments which are supposed to be real.
    2. Accepted synopses are z = feval(x,y,f) and z = feval(x,f) that support only one output argument. So,[Pmag, Pph_rad] = feval(P, polar); can't be a legal call.

    For your user case, since polar works only in a matricial way that is not the required one, just avoid it:

    P = 1 ./ (-10 * w.^2 + 0.1 * %i * w);
    [Pmag, Pph] = (abs(P), atan(imag(P),real(P)));
    

    About polar evolutions, do not hesitate to comment the reports #13486 and #16806.