|[eigen] uniform spaced cpx sinusoid via rootfinding, was Re: FFT module thoughts|
[ Thread 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:
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];
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)
fac = factors(k);
nth_root = y(2)^( 1/fac); % Use root-finding here
y_mult = [1 ];
raised = 1;
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)