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

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)
```

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