0

If I want to add 3dB per octave to signal (e.g. to flatten out the power spectrum of an exponential sine sweep), is it really just as simple as...

  1. Take the Fourier transform of the time-domain signal, let's call this A(f).
  2. Multiply this by the square root of the frequency (square root because power goes like amplitude squared)?
  3. Maybe normalize by the lowest frequency, so that there's unity gain for that frequency.
  4. Take the inverse Fourier transform to get back to the time domain.

Mathematically, the filter would be, simply,...

A'(f) := A(f) * sqrt ( f / f_low )

Is this right, and/or is there a better way?
(I've searched around and...it seems this is such a simple matter that people don't post how to do this.)

sh37211
  • 101
  • 7
  • See: http://dsp.stackexchange.com/questions/6220/why-is-it-a-bad-idea-to-filter-by-zeroing-out-fft-bins – MBaz Dec 04 '15 at 19:27
  • The post you refer to is about zeroing out bins. How is that relevant? I'm not talking about zeroing out bins. – sh37211 Dec 04 '15 at 19:41
  • you could also consider increasing the amplitude of your expontential sine sweep as the frequency increases. $$ $$ in fact, why do you even want to do this? if you're doing an exponential frequency sine sweep (as opposed to linear frequency sweep), ain't it because you want equal energy per octave instead of equal energy per Hz? – robert bristow-johnson Dec 04 '15 at 20:45
  • A whole international industry of people would want to do this because exponential sine sweep is the modern method for constructing acoustic impulse responses for measuring reverberation time (e.g., Farina, 2006 -- See http://pcfarina.eng.unipr.it/Public/Papers/226-AES122.pdf). This method is the international standard method, ISO 3382-2:2008. The exponential gives you a regular pattern of pre-echos which can be easily studied. The downside is that there is a -3dB/octave roll off from the exponential sweep which much be corrected for in order to render true ("flat") impulse results. – sh37211 Dec 04 '15 at 22:03
  • The signal is generated with uniform amplitude because of loudspeaker performance. The +3dB correction is intended as a post-processing correction effect. (Apart from that, the test signals have already been recorded and gaining access to the room again is not feasible. ) – sh37211 Dec 04 '15 at 22:07

2 Answers2

1

take any of the pinking filters you can find and swap the zeros with the poles.

e.g. here's my pinking filter with 3 real poles and 3 real zeros (all inside the unit circle):

$$ H(z) = A \frac{(z-q_1)(z-q_2)(z-q_3)}{(z-p_1)(z-p_2)(z-p_3)} $$

$p_1$ = 0.99572754, $q_1$ = 0.98443604

$p_2$ = 0.94790649, $q_2$ = 0.83392334

$p_3$ = 0.53567505, $q_3$ = 0.07568359

just swap the poles with the zeros.

robert bristow-johnson
  • 20,661
  • 4
  • 38
  • 76
  • I have used this filter before and found it and its inverse to to work well. From r b-j's notes elsewhere: "the response follows the ideal -3 dB/octave curve to within + or - 0.3 dB over a 10 octave range from 0.0009nyquist to 0.9nyquist" – Olli Niemitalo Dec 04 '15 at 23:46
  • yeah, it's cheap (inexpensive) and not particularly great. somewhere else i said, if i were to ever do it again, i would make it 5 poles, 4 zeros. never did it again. and i can't remember what value to put in for $A$. – robert bristow-johnson Dec 05 '15 at 02:19
  • I minimax optimized it for the same band for 5 poles, 4 zeros: $q_1=0.698258$, $p_1=0.378332$, $q_2=0.937174$, $p_2=0.862595$, $q_3=0.985792$, $p_3=0.970548$, $q_4=0.996652$, $p_4=0.993022$, $p_5=0.998655$. – Olli Niemitalo Mar 09 '16 at 10:47
  • (Poles and zeros were constrained to be real, dunno if that was too smart) – Olli Niemitalo Mar 09 '16 at 21:40
  • hey @OlliNiemitalo, you got 2668 rep. why don't you self-nominate yourself (that's evidently the only way, i can't nominate you) for moderator? – robert bristow-johnson Mar 10 '16 at 21:24
  • and thanks for doing that pinking filter. i will try it out. so you're saying that you minmaxed the deviation from $\sqrt{f}$ for the band from 0.0009 Nyquist to 0.9 Nyquist? if i can suggest something, i might suggest trading a little of the super-good accuracy in that 10-octave band for a little larger band, maybe 12 octaves. from 0.0002 Nyquist to 0.9 Nyquist. then the pinking filter will be good for up to 192 kHz sample rate. and i do think it's a good idea to constrain to real poles and zeros that are alternating. meanwhile i'll be copying your numbers to MATLAB. – robert bristow-johnson Mar 10 '16 at 21:30
  • Yes, minimax abs error in magnitude from 20/22050 Nyquist to 20000/22050 Nyquist to be precise. The max abs error with 5-pole 4-zero is 15 % compared to 3-pole 3-zero. I read that they will remove the moderator status from inactive moderators and I don't want to moderate regularly. I can leave the optimization overnight for a more wide-band filter. – Olli Niemitalo Mar 10 '16 at 21:40
  • if you redo this, why not from 20/96000 to 20000/22050 Nyquist? or maybe from 20/48000 to 20000/22050 Nyquist? so you're saying at present the deviation is about +/- 0.05 dB within the 10-octave band? – robert bristow-johnson Mar 10 '16 at 21:46
  • 20/96000 to 20000/22050 Nyquist sounds fine. I didn't look at the dB error nor didn't I optimize for it. – Olli Niemitalo Mar 10 '16 at 21:52
  • well, when the error is very small, whether it's dB or linear gain doesn't matter much. you have to find the $A$ factor so that the overall gain of your pinking filter lines up with the $\frac{1}{\sqrt{f}}$ function that you're matching it against. i'm just about to copy your numbers to MATLAB and check this out. – robert bristow-johnson Mar 10 '16 at 21:54
  • Well I have much more error in dB scale at high frequencies. Around 0.00037 dB at 20 Hz and 0.37 dB at 20 kHz. Maybe I should optimize the error in dB? – Olli Niemitalo Mar 10 '16 at 22:11
  • how're you optimizing? i thought you said minmax? oh, i get it. since the gain is 30 dB less at 20 kHz, then the error might be a lot more. ya, you should minimize the maximum dB error over the frequency range of interest. that means there will be a gain on your $\frac{1}{\sqrt{f}}$ curve that will place in right in the middle of your filter curve with equal-size + and - error. – robert bristow-johnson Mar 11 '16 at 01:32
  • I'm using home-brew differential evolution. Seems that for minimax ±0.1 dB error over 20/96000 Nyquist to 20k/22050 Nyquist the sweet spot is 4-real-pole 5-real-zero: $\bar q = {-0.2213191,0.4563575,0.8971910,0.9844598,0.9977857}$, $\bar p = {0.7506480,0.9596204,0.9940794,0.9992640}$ The coefficient $A=0.930487\sqrt{\omega}$ for unity gain at frequency $\omega$ in radians. – Olli Niemitalo Mar 12 '16 at 15:21
  • feel free to add that to the answer, after "firing up matlab" if you wish – Olli Niemitalo Mar 12 '16 at 19:53
  • whoa, you have a negative zero. we'll have to put another pole at z=0 to make this realizable. – robert bristow-johnson Mar 12 '16 at 20:16
0

If you are trying to convert the sweep response to an impulse response, you can do this directly by dividing the discrete Fourier transform (DFT) of the sweep response by the DFT of the sweep, and by taking the inverse DFT of the result. The division automatically compensates for the "pink" spectrum of the exponential sweep. Avoid division by zero and near-zero values at the very lowest and at the very highest frequencies; ramp the ends of the DFT of the sweep to magnitude of 1.

Olli Niemitalo
  • 13,491
  • 1
  • 33
  • 61