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

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

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
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];
    % inflate the vector by a factor of "fac" via outer product
    y = reshape( y_mult * y, 1,length(y)*fac );

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

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

Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/