Convert frequency response to impulse response

Question

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?)

Answer 1

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);

Answer 2

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.