FFT normalization factor depends on FFT frame size (which is equal to blocksize for fft~ & Co.), overlap, and window type.

In my experience you can normalize where ever you want: in time domain before FFT, in frequency domain, or in time domain after FFT. It is also possible to normalize in two steps, each of sqrt(normalization-factor), for example when you want the bin content normalized for certain processing or analysis jobs. Then you do one step before processing and the second after processing.

In frequency domain, you can either normalize amplitude of polar coefficients, or real and imaginary part of rectangular coefficients. This has no consequence for normalization factor.

If you would do normalization outside the reblocked subpatch, that is also possible, before or after FFT patch. But you must still consider blocksize and overlap.

You can verify all these possibilities by tweaking mentioned patch. The fact that normalization factor does not depend on the point in the process where you apply it, is convenient. However, for analysis or processing, that point of normalization does matter sometimes.

Consider for example a single-sample click: it's energy is spread over all frequency bins. With no normalization before spectral analysis, magnitude would be 1 for each bin. Next, consider a perfectly white noise, or a complex linear chirp over the full spectrum. With sqrt(normalization-factor) before analysis, the magnitude is 1 for each bin. Next, consider a pure sinusoid. With full normalization before analysis, the magnitude is 1 for a single bin, or the leaked equivalent of that.

You see, it depends on signal type how spectral magnitude can be shown independent of FFT size. Probably, the signals you want to process resemble a white noise rather than single pulses or full-scale pure sinusoids. So my guess is, that you would opt for the two-step normalization if you want to mix spectral data of different FFT sizes.

Good luck,

Katja

You actually should upsample, lowpass, distort, and lowpass again. The spectra of digital signals are periodic; it's technically not limited to the sample rate. The sample rate determines the size of the period. For real signals, you have frequencies from 0 to the Nyquist frequency, and then everything between Nyquist and the sample rate is a mirror image of the spectrum below Nyquist (you could think of them as the aliased frequencies). That defines one period, and it gets repeated further up the frequency range. In other words, you have the spectrum from 0 to SR, and that gets repeated at SR to 2*SR, and again at 2*SR to 3*SR, and so on.

Now, when you upsample, the parts of the spectrum above the original Nyquist will fall below the new Nyquist. And when you send it through [tanh~], those frequencies will produce new frequencies, some of which will alias in the new sample rate, and some of which will fall below Nyquist when you downsample back to the original sample rate. You probably don't want that. So, you'll need to filter after you upsample to remove the repeated spectrum. And the [tanh~] will produce so many partials that you'll need to filter again before you downsample.

I would recommend upsampling by a factor of at least 8. You'll still get some aliasing, but what gets aliased will probably get masked. I think Pd only lets you upsample by powers-of-two, so the next one would be a factor of 16. That should be high enough, though it could be cpu intensive. As for the filters, they should just be far enough below the original Nyquist so most or all of what is above it is filtered out. I would recommend using [lp10_cheb~] for this as it has a very steep roll-off. You could probably set it to about 18kHz without aliasing.

copy below dots in notepad etc...safe as ansi name.pd
........................................................................................................................
#N canvas 413 15 740 910 10;
#N canvas 559 52 558 609 fft 0;
#X obj 19 61 inlet~;
#X obj 195 217 inlet;
#X obj 29 92 rfft~;
#X obj 29 125 *~;
#X obj 60 125 *~;
#X obj 29 155 sqrt~;
#X obj 29 181 biquad~ 0 0 0 0 1;
#X text 93 93 Fourier series;
#X text 98 146 magnitude;
#X text 96 131 calculate;
#X text 21 3 This subpatch computes the spectrum of the incoming signal
with a (rectangular windowed) FFT. FFTs aren't properly introduced
until much later.;
#X text 83 61 signal to analyze;
#X text 193 164 delay two samples;
#X text 191 182 for better graphing;
#X obj 231 236 inlet;
#X text 284 234 toggle to graph repeatedly;
#X text 262 212 bang to graph once;
#X obj 19 295 tabwrite~ E09-signal;
#X obj 231 298 tabwrite~ E09-spectrum;
#X obj 29 205 /~ 4096;
#X msg 195 322 \; pd dsp 1;
#X obj 231 259 metro 70;
#X obj 332 109 block~ 4096 1;
#X obj 31 237 *~ 10;
#X connect 0 0 2 0;
#X connect 0 0 17 0;
#X connect 1 0 17 0;
#X connect 1 0 18 0;
#X connect 1 0 20 0;
#X connect 2 0 3 0;
#X connect 2 0 3 1;
#X connect 2 1 4 0;
#X connect 2 1 4 1;
#X connect 3 0 5 0;
#X connect 4 0 5 0;
#X connect 5 0 6 0;
#X connect 6 0 19 0;
#X connect 14 0 20 0;
#X connect 14 0 21 0;
#X connect 19 0 23 0;
#X connect 21 0 17 0;
#X connect 21 0 18 0;
#X connect 23 0 18 0;
#X restore 50 125 pd fft;
#X obj 109 68 bng 18 250 50 0 empty empty empty 0 -6 0 8 -262144 -1
-1;
#X obj 109 89 tgl 18 1 empty empty empty 0 -6 0 8 -262144 -1 -1 1 1
;
#N canvas 0 0 450 300 (subpatch) 0;
#X array E09-signal 882 float 0;
#X coords 0 1.02 881 -1.02 200 80 1;
#X restore 207 18 graph;
#N canvas 0 0 450 300 (subpatch) 0;
#X array E09-spectrum 259 float 0;
#X coords 0 0.51 258 -0.008 259 130 1;
#X restore 179 129 graph;
#X text 216 104 ---- 0.02 seconds ----;
#X text 271 147 SPECTRUM;
#X obj 24 512 gemwin;
#X obj 24 467 tgl 15 0 empty empty empty 17 7 0 10 -262144 -1 -1 1
1;
#X msg 50 469 create;
#X msg 53 489 destroy;
#X obj 191 288 gemhead;
#X obj 191 339 t a a a;
#X obj 180 368 GEMglEnd;
#X obj 236 367 GEMglBegin;
#X obj 362 350 GLdefine GL_LINES;
#X obj 362 320 loadbang;
#X obj 341 321 bng 15 250 50 0 empty empty empty 17 7 0 10 -262144
-1 -1;
#X obj 209 508 gemlist;
#X obj 209 429 f 256;
#X obj 209 448 until;
#X obj 208 472 t b b;
#X obj 293 470 + 1;
#X obj 210 396 t b b a;
#X obj 276 446 f 0;
#X obj 262 469 f;
#X obj 210 547 GEMglVertex2f;
#X obj 210 603 GEMglVertex2f;
#X obj 295 560 * -1;
#X obj 260 526 - 4;
#X obj 191 312 translateXYZ;
#X floatatom 242 288 5 0 0 0 - - -;
#X obj 295 520 tabread E09-spectrum;
#X obj 260 506 / 32;
#X obj 52 12 adc~;
#X obj 51 39 hip~ 5;
#X obj 51 65 *~ 1;
#X floatatom 94 42 5 0 0 0 - - -;
#X obj 210 579 spigot;
#X obj 149 564 tgl 15 0 empty empty empty 17 7 0 10 -262144 -1 -1 0
1;
#X text 404 -74 i think its a bouchard patch;
#X text 66 578 just for the symetry;
#X connect 1 0 0 1;
#X connect 2 0 0 2;
#X connect 8 0 7 0;
#X connect 9 0 7 0;
#X connect 10 0 7 0;
#X connect 11 0 30 0;
#X connect 12 0 13 0;
#X connect 12 1 23 0;
#X connect 12 2 14 0;
#X connect 15 0 14 1;
#X connect 16 0 15 0;
#X connect 17 0 15 0;
#X connect 18 0 26 0;
#X connect 19 0 20 0;
#X connect 20 0 21 0;
#X connect 21 0 18 0;
#X connect 21 1 25 0;
#X connect 22 0 25 1;
#X connect 23 0 19 0;
#X connect 23 1 24 0;
#X connect 23 2 18 1;
#X connect 24 0 25 1;
#X connect 25 0 22 0;
#X connect 25 0 33 0;
#X connect 25 0 32 0;
#X connect 26 0 38 0;
#X connect 28 0 27 2;
#X connect 29 0 26 1;
#X connect 29 0 27 1;
#X connect 30 0 12 0;
#X connect 31 0 30 1;
#X connect 32 0 26 2;
#X connect 32 0 28 0;
#X connect 33 0 29 0;
#X connect 34 0 35 0;
#X connect 35 0 36 0;
#X connect 36 0 0 0;
#X connect 37 0 36 1;
#X connect 38 0 27 0;
#X connect 39 0 38 1;

Patches IReditor and IRanalyser are added to the ExpoChirpToolbox.

The editor is a basic sample-editor for sample-precise trimming of the IR, windowing the tail(s), and normalising the maximum spectrum magnitude to 0 dBFS (full scale). The analyser presents spectral info of an edited IR, and you can test the IR as FIR filter with test signal, sound file or soundcard input. Later, the analyser will be extended to do parameter extraction, but since this will again require a lot of research, I want to share the patches now. Together with [expochirp-generator] and [IR-recorder], they already form a complete toolchain for generating your own FIR filters from acoustic environments.

Katja

ps .zip file attached this time (many thanks to the people who repaired this forum feature)

http://www.pdpatchrepo.info/hurleur/ExpoChirpToolbox.zip

@bassik said:

In this case I don't really get the point why you are researching the exact sample value of the peak value. In my experience, analysis of simple IR don't give substantial different results if the time zero of the IR is not exactly on the peak

Some more experimenting has learned me that it depends on the system under test. For an IR of the MacBook internal soundcard loopback, the problem is manifest. For an IR of my room, that is not the case. The difference is probably this: the soundcard's spectrum is near to perfect and there is only a direct response. The room test response missed a lot of the low spectrum end (also due to lousy 'test equipment'), and the direct response is buried in reflections because of the small room size. It's funny that the room response is much easier to trim than the loopback response.

Searching the internet I found the name for the phenomenon: fractional delay. And of course implementations for fractional delay, which can also be used to get rid of fractional delay. As far as I can see now, we'll only need it for cases where the direct response has a prominent role. For example, when we want to produce inverse filters for correction of spectral deficiencies in test equipment. We're not at that point yet.

Anyway, I did a patch with a model of fractional delay, which clearly illustrates the problem. (See attachment). And if you do a loopback test with [IRrecorder], you'll see what I mean as well.

edit: bah, again impossible to attach file. (Did I exceed my upload quotum?). It's here:

http://www.katjaas.nl/temp/fractional-delay.pd.zip

Katja

@katjav said:

Hey Bassik, thanks so much for your test results.

Apart from the 125 Hz the found RT are comparable indeed. But the magnitude plots of the second case? Should M0002_S01_R02.etm be comparable to Test 2b proc.wav? Even without a reference, the figure should show a similar shape because it is in deciBels, no?

But well, at least it shows that the Pd patch is not complete nonsense. We're on the right track.

I wonder how the magnitude of a full IR is computed in EASERA. For a 20 second chirp the IR it is at least 40 seconds, anticausal and causal part taken together. Or is the analysis limited to the causal part? No, I would think that full IR means the whole thing with the left side included. Pd can not do such long FFT's, 131072 points is the maximum. Probably the magnitude could be computed as an average of several overlapping frames. But the phase information is lost then.

I consider writing an FFT object for Pd that operates on a named buffer, and works with 64 bit floats internally. Maybe it is then possible to do larger FFT's. This is also important for calculating inverse filters (correction filters), since the full IR spectrum is sometimes needed for that. There is also something wrong with [fft~] and [rfft~] phase interpretation: it does not do a centered FFT. So with an IR centered you get these alternating phases instead of a zero-phase representation. Therefore, complex multiplication of spectra will give wrong results. Quite annoying. A new FFT object should have an option to select centered and non-centered FFT so you could work well with centered IR, and with causal part only. Hmm, plans....

Katja

Hello Katja,

Magnitude is expressed in dBFS and for definition our recorded signal should not go above 0 dBFS hence clipping in the analog domain.
So there is something wrong in EASERA when computing our IRs; I haven't figure out yet as the manual does not contain any reference to wav file loading and analysis.

Regarding which part is computed in Easera:
For acoustic measurements (RT, STI etc) only the casual part should be used.
The purpose of the swept sine technique is to get all the reproduction system non-linearities in the anti casual part (leftmost part) of the deconvoluted IR in order to then eliminate them.
So all the analysis should be done on a IR that has the time 0 where the peak of the IR is; also it will be then windowed to eliminate the late part that is not useful to be investigated as it might contain other artifacts of the deconvolution process due to e.g. the noise in the environment during measurements.

I hope I have understood your question.

Glad this keeps you thinking about new objects for PD.

I will review your post on the new chirp formula.

Thank you

Bassik

Bassik, there may be some confusion. This last patch chirp~test.pd is only a patch to test the qualities of the generated chirp itself.

I am kind of swallowed by details of the method. Instead of quickly producing a working IR measurement tool, I am now trying to refine the method, starting with refining the chirp itself. Last time I mentioned the sensitivity for parameter settings. Now I have found a systematic method to optimize a chirp. The formula is different from the regular ones as found in texts. This has lead to yet another chirp object [expochirp~] (yes it's exponential, not logarithmic, also see Farina's 2007 article).

The regular exponential chirp formula's make a chirp start at sine phase to avoid a step, but they do not care about the phase at the end of the function, therefore it may well end with a step and produce excessive ripple at the high spectrum end. This is traditionally relieved by windowing, or by manually cutting the chirp at the last zero-crossing as Farina suggests. My new formula computes the chirp in such a way that it will start and stop in sine phase. No need for windowing or cutting, and the ripples are minimized.

Still there will always remain ripples in a chirp's spectrum, and I am now trying to produce inversion-filters to clean up a system's measured IR. For the high spectrum end, this is fairly well possible. But in the low spectrum end, where it is more important, the ripple character of exponential chirps is such that conventional methods fail. I did find some solution but it's not without compromise. Therefore I am investigating alternative methods.

I am sorry that proceedings are slow now. But my hope is to contribute to a good quality tool in the end. I'll post details on the new chirp formula soon.

Katja

I would like to create an RGB color picker based on the the XY color spectrum. I have some LEDS and would like to select the color using the spectrum rather the three faders. this would require controlling the 3 faders using only 2.

I created a patch but could only get this far.

any help would be greatly appreciated. Thank you in advance.

here is a link of what i'm talking about.

http://www.google.com/imgres?imgurl=http://blommit.com/wp-content/uploads/2008/12/spectrum.jpg&imgrefurl=http://blommit.com/2008/12/27/the-ones-we-cant-see/&h=400&w=600&sz=21&tbnid=HguLqKMqh9ywVM:&tbnh=90&tbnw=135&prev=/images%3Fq%3Dcolor%2Bspectrum&zoom=1&q=color+spectrum&usg=__iomySBgMqPG4XqpLHr0-rpOReG8=&sa=X&ei=PP1ZTYqNNcTYgAe4zLCdDQ&ved=0CCYQ9QEwAA

http://www.pdpatchrepo.info/hurleur/X_Y_spectrum_picker.pd

Stupid me. I overlooked the fact that magnitude inversion has to be done in one single fft frame as well. For your test signal of 352800 samples that means a 2^19 point fft. Though it need not theoretically be a problem, Pd can not handle fft with that blocksize.

Several other options remain.

1. If you could do with a linear sweep test signal, you don't need to invert magnitude spectrum because a carefully constructed linear sweep has flat spectrum. The time-reversal to invert phases would suffice.

2. Invert the log sweep's spectrum in another software that can do very long fft's.

Will [partconv~] be able to handle long impulse responses? According to the object's c code it can do up to 256 partitions. Your test signal is done in 178 partitions of size 2048 so that should be ok.

What is by the way the formula of your test signal? Did you create it in a Pd patch? I'd like to know.

Katja

@lead said:

maybe process the lists from cartesian to polar coordinates (or visa versa) then *-1 +1?

Not sure what you mean exactly. Admittedly, my suggestion was rather cryptic as well.

Anyway, I was thinking about the situation of a test signal: how long is that signal? Probably too long to fit in one fft frame of reasonable length. In that case, I don't think you can invert it's phase spectrum in frequency domain. But no problem: just do a time reversal in signal domain, and it's phase spectrum will be inverted automatically. To be explicit here: I mean you must time-reverse the samples of the original test sweep (not the test result signals). Time reversal does not change the magnitude spectrum, only the phase spectrum.

So remains inverting the magnitude spectrum of the time-reversed test sweep. A log sweep does not have a flat magnitude spectrum, so you need to compute the multiplicative inverse of the magnitude coefficients (in frequency domain, using polar coordinates indeed here). Beware of possible zero-magnitude points in the original, they can not be inverted of course. After magnitude inversion, go back to cartesian coordinates and from them to time domain. Now you have your deconvolution kernel in time domain.

The deconvolution kernel is the test sweep with phase and magnitude spectrum inverted. It can be considered a (very long) FIR filter. With this kernel you can do fast convolution of the test results, to deconvolve these. For this, you can use Ben Saylor's [partconv~] object. It does 'partitioned convolution' in frequency domain, with zero padding to avoid circular convolution.

Good luck. Please let us know if you have succes. It is an interesting technique to create impulse responses of acoustic conditions. Expensive reverb simulators probably rely on similar techniques. What is actually your purpose for it?

Katja