1
$\begingroup$

I am running an acoustic simulation. I put in a frequency and get out a complex number; the frequency response for that frequency.

I wish to generate a 1024-sample impulse response (supposing my output sample-rate is 44.1kHz)

How do I go about it?

I think I have the first step: take the Nyquist frequency 22.05kHz, and cut it into 512 equal slices. So each slice will have a bandwidth of 22050/512 Hz

Now I feed these frequencies into the simulator and get back 512 complex numbers; the first will be 0+0i corresponding to 0Hz

Now I'm following an instruction that says: 'consider 0Hz now to be the symmetry point, and add 511 negative frequencies.

So I have { Re(511) + i.Im(511), ..., Re(1) + i.Im(1), 0 + 0i, Re(1) + i.Im(1), ..., Re(511) + i.Im(511), }

1023 frequencies.

Now I'm confused. I don't know what I'm doing any more.

Apparently I need to IFFT this, and that is going to give me an 1023 sample impulse response?

Could somebody help me get my feet back on the ground?

What would the MatLab code for this look like?

(Also, is there any way to do it from the Bash shell using Python or something similar?)

$\endgroup$
  • $\begingroup$ Yes, you can use an IFFT on the conjugate symmetric vector (You might want to zero-pad it on both sides first). Why is this confusing? $\endgroup$ – hotpaw2 Jan 19 '14 at 17:23
  • $\begingroup$ If it was a steady signal, and I know the frequency components, then I can reconstruct the signal in time domain by superposition. I can think of each component as a sinusoid, a point moving round a circle on the origin of the complex plane. But I can't get from here to that result. Do I have to abandon this line of thinking? $\endgroup$ – P i Jan 19 '14 at 18:43
  • $\begingroup$ An IFFT is a fast exact (within numerical rounding) equivalent computation to constructing by each sinusoidal component. $\endgroup$ – hotpaw2 Jan 19 '14 at 20:46
3
$\begingroup$

Here is a MATLAB example. Inefficient code but hopefully easy to read

%% Example: first order bandpass from 200 to 2000 Hz sampled at 44.1 kHz
n = 1024; % FFT length
fs = 44100; % sample rate
f0 = [200 2000]; % bandpass range
[b,a] = butter(1,2*f0/fs); % these are the difference equation coefficients
b0 = b(1); b1 = b(2); b2 = b(3);
a0 = a(1); a1 = a(2); a2 = a(3);  % for easier readability

% Step 1: build normalized frequency vector
om = 2*pi*(0:n/2)'/n;

% step 2: calculate the frequency response at every frequency
z = exp(j*om); % this is the "z" of the z-transform
% calculate the Z transform , 
H = (b0+ b1*z.^-1 + b2*z.^-2)./(a0+ a1*z.^-1 + a2*z.^-2);

% Step 3: make it conjugate symmetric
H = [H; conj(H(end-1:-1:2))];

% Step 4: inverse FFT
h = real(ifft(H));

% Step 5: verify against direct evaluation of the time domain difference equation
d0 = zeros(n,1); d0(1) = 1; % unit impulse
href = filter(b,a,d0);
err = 10*log10( sum((h-href).^2)./sum(href.^2));
fprintf('Error = %6.2f dB\n',err);
$\endgroup$
1
$\begingroup$

I've put up an answer at https://stackoverflow.com/questions/21209017/python-convert-frequency-response-to-impulse-response

Basically if my frequencies are:

{ f_dc, f_1, ..., f_k, f_nyquist }  

    (f_dc=0, f_nyquist for me = 22050 a.k.a freq limit for 44.1kHz audio)

I have to construct:

{ DC, f_1, ..., f_k, nyquist, f_k*, ..., f_1* }

    with DC = nyquist = 0 (note: z* is the complex conjugate of z)

Then I iFFT.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.