[eigen] uniform spaced cpx sinusoid via rootfinding, was Re: FFT module thoughts

[Mark Borgerding <mark@xxxxxxxxxxxxxx>]

To create a root-finder, I'd probably start by trying to use 
http://en.wikipedia.org/wiki/Newton%27s_method.  There might be a more 
suitable algorithm for finding complex roots of unity.

One could use this root-finder to create an entire cycle via the 
algorithm in the attached octave/matlab script.
The idea is similar to http://en.wikipedia.org/wiki/CORDIC

The only reason I can think of not to do this for the general case is 
fear of worst-case numerical performance (e.g. n is a large prime).  One 
thought: This  might be mitigated by comparing the complex sinusoid with 
its time-reversed and conjugated self.


On 04/05/2012 11:02 AM, Pavel Holoborodko wrote:
> Hi Mark.
[...]
> I would really appreciate if you would provide some links on how to re-write make_twiddles using root finding approach.

clear all;close all;

factors=[2 2 3 5];
n=prod(factors);% 60
theta=(2*pi*(0:(n-1))/n)';
y_desired = exp(j*theta);

figure;plot(theta,[real(y_desired) imag(y_desired)] );legend('real','imag')

% TODO: don't presume factors(1)==2, i.e. n is even
y = [1 -1];% cexp(j0) cexp(jpi)

for k=2:length(factors),
    fac = factors(k);
    nth_root = y(2)^( 1/fac); % Use root-finding here
    y_mult = [1 ];
    raised = 1;
    for m=2:fac,
        raised = raised * nth_root;
        y_mult = [y_mult;raised];
    end
    % inflate the vector by a factor of "fac" via outer product
    y = reshape( y_mult * y, 1,length(y)*fac );
end
y=y(:);

figure;plot(theta,[real(y) imag(y)] );legend('real','imag')

rmse = norm(y_desired-y)/norm(y)