Setting fft bin real and imaginary based magnitude and phase

DSP, Plug-in and Host development discussion.
RELATED
PRODUCTS
User avatar
SNFK
KVRist
Topic Starter
70 posts since 8 Nov, 2020

Post Wed Jan 19, 2022 6:14 pm

I know that magnitude = sqrt(real*real+imaginary*imaginary) and phase = -atan(imaginary, real). I want to be able to set the FFT bin based on either a new magnitude or phase. I ended up with this function to get the reals, but it doesn't seem to work:

re = sqrt(mag^2-im^2)

Any suggestions?

User avatar
aciddose
KVRAF
12500 posts since 7 Dec, 2004

Post Wed Jan 19, 2022 6:33 pm

As far as I understand it the decision about the imaginary phase angle is arbitrary. You need to remember a few key points about discrete Fourier transforms:
  • The transform only applies to continuous (infinite length) cyclic signals (like sin())
  • There is absolutely no difference between a discrete Fourier transform and additive synthesis
Given those two points the obvious solution should occur to you: transform a cycle of an arbitrary waveform and examine the complex output. If the output contains any imaginary phase component it makes the process of synthesizing the complex table for inverse transform non-trivial and in most cases fractal (seems arbitrary.)

Where you have very simple additive tables as input (sine, triangle, pulse, ramp) you should observe pure real output and see no difference between additive and the inverse Fourier transform. (The real output is 'rotated' into the imaginary frequency value. All 'real' frequency values should approximate zero for a 'zero phase' input.)

The complex output from a Fourier transform of a real input is really just an additive pattern packed in a particular way. Given an input of 65536 length with a single cycle of sin(), the imaginary frequency value IF[1] should equal exactly length/2 or 32768. So the normalization factor is 2/length.

When you transform a continuous cyclic function with a rectangular (no) window you will get an output of exactly the amplitude of each harmonic partial in the input.

Maybe this might help:

Code: Select all

T amplitude[len/2] = { ... };
for (i = 1; i < len/2; i++) {
	RF[i] = T(0);
	IF[i] = amplitude[i - 1] * len/2;
	// mirror top half of complex values
	RF[len-i] = T(0);
	IF[len-i] = amplitude[i - 1] * -len/2;
}
That should work for 'zero phase'. I've never seen the point of using Fourier to do additive (I'm sure has optimization applications) so I'm not sure how phase works. I assume you just 2d rotate the amplitude into the complex according to the phase. So the phase value should actually just be a rotation for a [cos(phase), -sin(phase)] pair. Think of your amplitude being a normalized "up" vector and RF[] and IF[] are like v.x and v.y you're rotating the amplitude vector into 2d "complex" space.

Right? So amplitude/magnitude might be a real, but that's just a complex vector with a zero imaginary value. So if amplitude = 1, it's actually {1, 0} which is just a 2d vector directly "up".

In the mirrored half of the frequency values we're pointing "down", which is the same as a 180 degree (?) complex rotation. The thing I don't know is whether you rotate the mirrored values in the same or opposite direction.

Image

Basically, calling it real/imaginary/complex is a load of bullocks and gibberish. It's not complex, it's very simple and all we're doing is rotating in 2d.
Free plug-ins for Windows, MacOS and Linux. Xhip Synthesizer v8.0 and Xhip Effects Bundle v6.7.
The coder's credo: We believe our work is neither clever nor difficult; it is done because we thought it would be easy.
Work less; get more done.

hugoderwolf
KVRist
227 posts since 1 Apr, 2009 from Hannover, Germany

Post Thu Jan 20, 2022 12:42 am

re = mag * cos(phase)
im = mag * sin(phase)

User avatar
Music Engineer
KVRAF
4079 posts since 8 Mar, 2004 from Berlin, Germany

Post Thu Jan 20, 2022 1:44 am

right. and also:

phase = atan2(im, re)

...you've got that formula wrong, too. mag = sqrt(re^2 + im^2) is correct

User avatar
aciddose
KVRAF
12500 posts since 7 Dec, 2004

Post Thu Jan 20, 2022 3:26 pm

hugoderwolf wrote: Thu Jan 20, 2022 12:42 am re = mag * cos(phase)
im = mag * sin(phase)
That clearly doesn't work on its own because it's not taking into account the mirroring of the upper half.

Are you saying this is correct?:

Code: Select all

// linear magnitude
const T magnitude[len/2] = { ... };
// phase in units of pi "radians"
const T phase[len/2] = { ... };
const T norm = T(len/2); // normalization_factor
for (i = 1; i < len/2; i++) {
	const T mag = magnitude[i - 1];
	const T phs = phase[i - 1];
	RF[i] = mag * cos(phs) * norm;
	IF[i] = mag * sin(phs) * norm;
	// mirror top half of complex values as complex conjugates
	RF[len-i] = RF[i];
	IF[len-i] = -IF[i];
}
Edited code to remove sign inversion on RF upper half as mentioned by mystran below.
Last edited by aciddose on Thu Jan 20, 2022 4:05 pm, edited 1 time in total.
Free plug-ins for Windows, MacOS and Linux. Xhip Synthesizer v8.0 and Xhip Effects Bundle v6.7.
The coder's credo: We believe our work is neither clever nor difficult; it is done because we thought it would be easy.
Work less; get more done.

mystran
KVRAF
6935 posts since 12 Feb, 2006 from Helsinki, Finland

Post Thu Jan 20, 2022 3:44 pm

aciddose wrote: Thu Jan 20, 2022 3:26 pm
hugoderwolf wrote: Thu Jan 20, 2022 12:42 am re = mag * cos(phase)
im = mag * sin(phase)
That clearly doesn't work on its own because it's not taking into account the mirroring of the upper half.
If you intend the time-domain signal to be real-only, then either (1) use a real-only FFT, which saves you about half the FFT computation or (2) set the "negative" frequencies to be the complex conjugates of the positive frequencies (ie. negate the imaginary part, but not the real part).
Preferred pronouns would be "it/it" because according to this country, I'm a piece of human trash.

User avatar
aciddose
KVRAF
12500 posts since 7 Dec, 2004

Post Thu Jan 20, 2022 7:13 pm

aciddose wrote: Wed Jan 19, 2022 6:33 pm As far as I understand it the decision about the imaginary phase angle is arbitrary.
Quoting myself to clarify here. I mean many things by this statement. In part I'm referring to the "handedness" of a system such as complex numbers, positive/negative or any N dimensional system like x,y or x,y,z.

Here's a quick google search result:
https://physics.stackexchange.com/quest ... transforms

What I've always found is the same answer: "Just because."

Basically, we just arbitrarily choose whether X and Y are width, height, depth in a system of two dimensions. Likewise whether a positive direction on an axis corresponds to leftward or rightward or upward or downward or forward or backward.

The choice about the imaginary phase angle of the real input signal is also completely arbitrary. We just choose [magnitude, 0], and then from there we choose that after the transformation we end with [0, magnitude], and that phase rotates counter-clockwise and so on.

It just seems weird that we had to choose at all. The same issue with a quadratic equation has two solutions/roots. And from there the issues of defining hyperspace, the uncertainty principle and so on. There are two solutions and we have to choose one of the two, we can't have both and there is no solution to these equations "in between".

Ultimately though what makes it seem confusing and arbitrary is the issue in itself. We are choosing arbitrarily and there is no higher system within which we can reduce this choice to only one option. This can be frustrating if you begin in search of some higher meaning through mathematics.

Here's a favorite quote of mine:
Bertrand Russell wrote: Pure mathematics consists entirely of assertions to the effect that, if such and such a proposition is true of anything, then such and such another proposition is true of that thing. It is essential not to discuss whether the first proposition is really true, and not to mention what the anything is, of which it is supposed to be true. Both these points would belong to applied mathematics.

We start, in pure mathematics, from certain rules of inference, by which we can infer that if one proposition is true, then so is some other proposition. These rules of inference constitute the major part of the principles of formal logic. We then take any hypothesis that seems amusing, and deduce its consequences. If our hypothesis is about anything, and not about some one or more particular things, then our deductions constitute mathematics.

Thus mathematics may be defined as the subject in which we never know what we are talking about, nor whether what we are saying is true. People who have been puzzled by the beginnings of mathematics will, I hope, find comfort in this definition, and will probably agree that it is accurate.
https://en.wikipedia.org/wiki/I_know_th ... ow_nothing
? wrote:The more I learn, the less I know.
Free plug-ins for Windows, MacOS and Linux. Xhip Synthesizer v8.0 and Xhip Effects Bundle v6.7.
The coder's credo: We believe our work is neither clever nor difficult; it is done because we thought it would be easy.
Work less; get more done.

hugoderwolf
KVRist
227 posts since 1 Apr, 2009 from Hannover, Germany

Post Thu Jan 20, 2022 11:50 pm

aciddose wrote: Thu Jan 20, 2022 3:26 pm
hugoderwolf wrote: Thu Jan 20, 2022 12:42 am re = mag * cos(phase)
im = mag * sin(phase)
That clearly doesn't work on its own because it's not taking into account the mirroring of the upper half.
My math is so good it would even work with FFTs of complex signals. 8)

Of course we can go on trying to answer questions that have never been asked in this thread, based on random assumptions on what the OP is actually trying to achieve (which we don't know). But the above answers the question that was asked.

DaveClark
KVRist
327 posts since 8 May, 2007

Post Fri Jan 21, 2022 11:39 am

aciddose wrote: Thu Jan 20, 2022 3:26 pm
That clearly doesn't work on its own because it's not taking into account the mirroring of the upper half.

Are you saying this is correct?:

Code: Select all

// linear magnitude
const T magnitude[len/2] = { ... };
// phase in units of pi "radians"
const T phase[len/2] = { ... };
const T norm = T(len/2); // normalization_factor
for (i = 1; i < len/2; i++) {
	const T mag = magnitude[i - 1];
	const T phs = phase[i - 1];
	RF[i] = mag * cos(phs) * norm;
	IF[i] = mag * sin(phs) * norm;
	// mirror top half of complex values as complex conjugates
	RF[len-i] = RF[i];
	IF[len-i] = -IF[i];
}
F(N-n)* = F(n) imposes restrictions on the full sequence of values for phase: setting the first half of the values for phase determines the other half of the values. If the F(n) are actually from a real-valued signal, then yes, the code or something very similar to it will work. Be extra careful with index calculations. See for example, William H. Press, et al. Numerical Recipes in C (FORTRAN, etc.): The Art of Scientific Computing. Second Edition, 1992, Cambridge University Press, pages 510 - 514. Their code is so-so, but their explanations are great.

User avatar
aciddose
KVRAF
12500 posts since 7 Dec, 2004

Post Fri Jan 21, 2022 2:30 pm

These threads often show up in search results for example so it doesn't hurt to lecture a bit. I hate finding a result that covers a particular question that is tangentially related to what I'm looking for but then doesn't go on to explain why or anything else about the area.

It's even more frustrating to try to identify these sorts of answers "why won't my mag/phase work with my fft?" = "you're using it as if it were a real-only inverse transform" = "what is a real-only fast inverse Fourier transform?", "what do I do if the library I use doesn't have this?"

Generally the information I'm looking for is: "How do I write this library myself from absolute scratch?", so having a bunch of esoteric mathematics jargon flung at me in a lecture is frustrating. As Carl Friedrich Gauss said in the quote I posted above, if it weren't for this jargon the underlying mathematics would be intuitive.

Understanding that the equations implicitly lead to the "imaginary" output storing the amplitude of partials from a zero-phase input signal in the "real" input because of an arbitrary "handedness" choice made hundreds of years ago is not at all intuitive! You'll also rarely find it stated so succinctly.

While the OP did not ask specifically it is implied by the format of the question asked that the OP also did not understand a lecture on discrete Fourier transforms encompassing the full suite of esoteric jargon.

Always remember these series of studies on related phenomena:
https://en.wikipedia.org/wiki/Dunning%E ... ger_effect
Wikipedia wrote:Moreover, competent students tended to underestimate their own competence, because they erroneously presumed that tasks easy for them to perform were also easy for other people to perform.
In other words experts in a field will severely underestimate the effort required to develop their expertise and have a tendency to assume others can just as easily comprehend the same fields. Information packed away in the "sub conscious" of the expert is neither present in the mind of the newcomer nor necessarily other experts.
Free plug-ins for Windows, MacOS and Linux. Xhip Synthesizer v8.0 and Xhip Effects Bundle v6.7.
The coder's credo: We believe our work is neither clever nor difficult; it is done because we thought it would be easy.
Work less; get more done.

mystran
KVRAF
6935 posts since 12 Feb, 2006 from Helsinki, Finland

Post Fri Jan 21, 2022 3:20 pm

aciddose wrote: Fri Jan 21, 2022 2:30 pm Understanding that the equations implicitly lead to the "imaginary" output storing the amplitude of partials from a zero-phase input signal in the "real" input because of an arbitrary "handedness" choice made hundreds of years ago is not at all intuitive! You'll also rarely find it stated so succinctly.
The zero-phase result is found in the real component (=cosines). It's just that if you're trying to synthesize waveforms then choosing the zero-phase gives you a waveform with one big spike at t=0 and instead you typically want to add together sines (=90 degrees out of phase), because this gives you a more sensible crest-factor (precisely because not all the partials are in phase with each other).

edit: also the sign of the imaginary component makes very little difference, because if you take the "forward" FFT and keep applying it then you just cycle through time-domain -> spectrum -> time-reversed-time-domain -> conjugate spectrum -> original-time-domain again...

edit2: for all intents and purposes, the FFT itself is a 90 degree rotation on a 2d plane where one axis is time and the other axis is frequency (and yes, this theoretically does generalize to "fractional fourier transforms" as well, though the actual math for those is a little bit too much for me).
Preferred pronouns would be "it/it" because according to this country, I'm a piece of human trash.

Return to “DSP and Plug-in Development”