Rather than use the default implementation of DFT I want to use a generalised version. Now obviously I could easily code it up myself using the formulas in those links but what I'm doing is horrifically computer intensive and my code is rarely 'efficient'. Nothing I code up will even approach the kind of efficient speed found in things like the implementation of DFTs found in the standard libraries of C, R, Mathematica, Matlab etc, as there are algorithms written by people who devote their entire research lives to doing numerical FTs quickly. Does anyone know which, if any, programming languages have a standard library for the generalised DFT? Preferrably I'd prefer to work with either R or Octave but beggars can't be choosers. I could do the middle ground of preprocessing a data block by multiplication of the entries by the appropriate roots of unity and then doing a standard DFT but that's still pretty inefficient. That seems to be the way the fastest Fourier transform in the West does it which I find a little unsatisfactory but then they are the FFTW! Hopefully someone here has some experience with them and can point me in the right direction. Thanks Please Register or Log in to view the hidden image!
I have no idea but I just wanted to make an uneducated guess at some possible links. Beggars can't be choosers...Please Register or Log in to view the hidden image! http://maltman.hmdc.harvard.edu/numal/resources/ http://en.wikipedia.org/wiki/List_of_numerical_libraries The Cooley-Tukey FFT Algorithm for General Factorizations
It doesn't look to me like the GDFT uses roots of unity, only - a and b aren't necessarily integers... Also the suggested preprocessing is just linear time (a single pass, in fact) in the input data, and that will be completely dominated by the O(n log n) for the FFT. If you want to do what the FFTW guys are talking about, I strongly suggest following their lead. (I originally wanted to suggest simply replacing the precomputed twiddle factors in a simple (iterative) FFT with the generalized ones, but I can't remember whether the FFT algorithm actually relies on property of the twiddles being roots of unity - it's been 5 or 6 years since I last looked at this stuff in detail...)
The data is something I construct so I'm able to define the size of the array to make a and b integers. The expression given on the Wiki page doesn't match that from the FFTW page, where they say just multiply an entry by \((-1)^{i+j+\ldots}\), unless I'm missing something?. I've just spent the last 2 hours or so messing around with Mathematica routines to try to get the right factors of \(e^{\frac{2\pi i}{N}\ldots}\) in an example array but there's a mistake somewhere as I get almost the answer I know it should be but it looks like its out by a single factor of something but its too late in the evening for me to be bothered now. I can see its basically doing a translation in each space, you shift the data, do a standard DST and then shift back and overall that combines into the expression on the Wiki page. Unfortunately I need sleep now.... While Googling I found a pdf talk on the DFT methods and about 1/2 a page in I thought "This seems familiar" before then realising the author had simply converted the Wiki page into a pdf! New backdrop but same headings, text, equations, pictures, even links! Talk about dishonest and lazy!
This looks relevant, but I don't have access to the article. http://dx.doi.org/10.1016/S0165-1684(97)00234-X Guoan Bi and Yanqiu Chen "Fast generalized DFT and DHT algorithms" Signal Processing 65 3, (March 31, 1998), Pages 383-390 Also likely hits: http://dx.doi.org/10.1016/0012-365X(85)90016-0 Thomas Beth "Generalizing the discrete fourier transform" Discrete Mathematics 56 2-3 (October 1985) Pages 95-100 http://dx.doi.org/10.1016/S0165-1684(99)00088-2 Vladimir Britanak and K. R. Rao "The fast generalized discrete Fourier transforms: A unified approach to the discrete sinusoidal transforms computation" Signal Processing 79 2, (December 1999), Pages 135-150 Related thesis on Fast FFTs. http://spiral.ece.cmu.edu:8080/pub-spiral/pubfile/thesis_112.pdf // Edit: These matlab people have made me think that this is easy. Is this much like what you were looking for? http://www.ee.columbia.edu/~marios/symmetry/sym.html Code: function X = gdft(x,a,b) % GDFT Generalized discrete fourier transform % X = gdft(x,a,b) % % The case of a=b=0 corresponds to dft % The case of a=1/2, b=0 corresponds to ofdft (aka odft) % % x: input signal, follows fft interface of frame per column % a: shift in frequency (usually 0 or 1/2) % b: shift in time (usually 0 or 1/2) % X: gdft of x % % Notes: This is the normalized gdft unlike the fft-ifft definition % in Matlab. This means that Parsevals theorem doesn't need % division by N % ------- gdft.m ------------------------------------------- % Marios Athineos, marios@ee.columbia.edu % http://www.ee.columbia.edu/~marios/ % Copyright (c) 2004-2005 by Columbia University. % All rights reserved. % ---------------------------------------------------------- % The default is the Odd Frequency DFT if nargin < 2; a = 1/2; end if nargin < 3; b = 0; end % Argument check if a<0 | a>1; error('Shift in frequency must be between 0 and 1 (inclusive)'); end if b<0 | b>1; error('Shift in time must be between 0 and 1 (inclusive)'); end % Get size of input signal [N,fnum] = size(x); % Pre weights (in time) nw = exp(-j*2*pi*a*[0:N-1]/N); nw = diag(sparse(nw)); % Post weights (in frequency) kw = exp(-j*2*pi*([0:N-1]+a)*b/N); kw = diag(sparse(kw)); % So the gdft is simply (normalized) X = full(kw*fft(full(nw*x)))/sqrt(N);
I know someone that actually has done some DFT's in C++, if you can use C++ let me know I'll see if he has ever finished that up. I think he built his own kind of like you are trying to do.
Turns out there's a 1 line code in Mathematica which implements precisely what I wanted. Rather than shifting the powers of unity by multiplying by fiddly factors you can just say RotateLeft[data.array,{a1,a2,...,aN}] and it'll cyclically permute the first index by a1 terms, the second by a2 etc. It actually moves the data in the array precisely like a translation and then after the DFT is done you just apply RotateLeft[DFT'd.data.array,{-a1,-a2,...,-aN}] and you're done Please Register or Log in to view the hidden image! Much better than the god awful code I was half way through writing! Rpenner that looks like the kind of thing I was thinking of, they've used the nice index manipulating features in Matlab to make it succinct. Mathematica doesn't have the same raw syntax power but fortunately there's plenty of nice protocols programmed in Please Register or Log in to view the hidden image! Thanks guys.
If you're looking to use Fourier transforms in a language which isn't solely dedicated to mathematical analysis... Python has either: SciPy - http://www.scipy.org/ (Or specifically the package within it...) Numpy - http://numpy.scipy.org/
I've got another question related to the numerical implementation of Fourier transforms and seeing as you guys have much more experience than I in them I thought I'd ask it here. My previous question related to shifting the data within the array/table its contained within. This was more for conceptual ease than mathematics. Given some array like Data[[a,b,c]] (I'm using Mathematica notation here but you're all good enough programmers to know what I'm talking about) I can take a 3D Fourier by just saying Fourier[Data]. Fine. I can take a 2D Fourier by specifying the array components, like Fourier[Data[[1]]] will do the 2D Fourier of the 2d array defined by Data[[1]]. I can take other 2D Fourier's which are parallel to the 'axes' by just doing something like Fourier[Data[[,1,]]] (okay, that's mixing Mathematica and R notation but you get my point). Now here's the kicker. What if I don't want to take a Fourier transform of a 'nice' 2D array but one which say goes diagonally across my 'cube' of information denoted by 'Data' ? It basically amounts to me having to rotate the data or pick out a rotated section of 'Data'. The problem is that when you do this you get rounding errors because I can only pick out integer elements, in the sense that I can't pick out the 2.5'th element of a list, only the 2nd or the 3rd. As a result when I then inverse Fourier transform back I don't what I wanted, I get some very crude approximation to it which is not fit for purpose. Are there methods which lessen this rounding error problem? Surely people have wanted to do Fourier transforms of rotated things before? Increasing the discretisation level of my data (ie taking loads more discrete samples from a continuous input) doesn't help because no matter your 'resolution' doing a line at 45 degrees through a square of pixels is a pain in the arse!! I need someone build a computer which works with continuous variables, not discrete approximations!! Please Register or Log in to view the hidden image!
Please Register or Log in to view the hidden image! AKAT-1, 1959 Yeah baby! See also Hybrid computer (Wikipedia): Hybrid computers are computers that exhibit features of analog computers and digital computers. The digital component normally serves as the controller and provides logical operations, while the analog component normally serves as a solver of differential equations. In general, analog computers are extraordinarily fast, since they can solve most complex equations at the rate at which a signal traverses the circuit, which is generally an appreciable fraction of the speed of light. On the other hand, the precision of analog computers is not good; they are limited to three, or at most, four digits of precision. Digital computers can be built to take the solution of equations to almost unlimited precision, but quite slowly compared to analog computers. Generally, complex equations are approximated using iterative numerical methods which take huge numbers of iterations, depending on how good the initial "guess" at the final value is and how much precision is desired. (This initial guess is known as the numerical seed for the iterative process.) For many real-time operations, the speed of such digital calculations is too slow to be of much use (e.g., for very high frequency phased array radars or for weather calculations), but the precision of an analog computer is insufficient. Hybrid computers can be used to obtain a very good but relatively imprecise 'seed' value, using an analog computer front-end, which is then fed into a digital computer iterative process to achieve the final desired degree of precision. With a three or four digit, highly accurate numerical seed, the total digital computation time necessary to reach the desired precision is dramatically reduced, since many fewer iterations are required. Probably not relevant, but I'm avoiding exam study.
What's this "cube" ? Confused about the data. Rotating a function rotates it's 2D Fourier transform. Regardless about how you rotate (in frequency or space domains) you wind up with pieces beyond the boundaries of the original square, which you can avoid by padding. If you are using typical convolution methods, this actually causes aliasing. http://morse.cs.byu.edu/450/lectures/lect18/ft-2d.slides.printing.6.pdf The following paper factors the shift and rotation of an inverse fft into : Z* .* conv( S .* f .* Z , Z* ) = Z* .* ifft( fft( S .* f .* Z) .* fft(Z*) ) where .* is element-by-element multiplication. http://afni.nimh.nih.gov/sscc/staff/rwcox/23dw.ps Other rotation ideas: http://www.cs.dartmouth.edu/reports/TR96-301.pdf \(\left( \begin{array}{cc}1 & -\tan \frac{\theta}{2} \\ 0 & 1 \end{array} \right) \left( \begin{array}{cc}1 & 0 \\ \sin \theta & 1 \end{array} \right) \left( \begin{array}{cc}1 & -\tan \frac{\theta}{2} \\ 0 & 1 \end{array} \right) = \left( \begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array} \right)\)
