How can I interpolate between two FFTs?

Question

I have two FFTs that represent two impulse responses. I want to interpolate between the two of them to generate a 3rd FFT in the manner shown in the diagram below. This is to be used in real time audio rendering filtering between two FFT filters.

Linear interpolation between each of the FFT's complex coefficients doesn't work as wanted.

FFT plots of two impulse responses 
(x = frequency/coefficient, y = frequency magnitude)

+                                   +                           
|                                   |                           
|        +                          |                           
|        |                          |                           
|        |                          |                           
|        |                          |                           
|        |                          |                           
|        |                          |                   +       
|        |                          |                   |       
|        |                          |                   |       
+--------+------------------+       +-------------------+-------+    
           Input A                              Input B

FFT plots of interpolation between A and B (Linear, and Wanted)
(x = frequency/coefficient, y = frequency magnitude)

+                                   +                                
|                                   |                                
|                                   |                                
|                                   |                                
|                                   |                                
|                                   |             +                  
|        +                          |             |                  
|        |                          |             |                  
|        |          +               |             |                  
|        |          |               |             |                  
+--------+----------+-------+       +-------------+-------------+    
Linear Interpolation (A+B/2)            Wanted Interpolation

How do I do the interpolation as described in the Wanted Interpolation graph?

Answer 1

This is the solution for a linear interpolation between two frequencies as a "weighted average". You would just need to decide what your weighting actually is in terms of magnitude or dB scale for the vertical axis. X is the interpolated frequency, weighted toward the stronger signal, and Y is the linear interpolation at the new frequency location of the signal strenths a and b.

Here is a simple example demonstrating the results that would be achieved with the formulas given:

Answer 2

Have a look at the script below. It does not solve all your problems, but it can be a starting point:

import scipy.interpolate

# 1) Here, I define some arbitrary initial spectrum, 
#    for example it consts of two paeaks
#    Here, I simplify and set the frequencies to be between 0 and 1 in 100 steps.
f = np.linspace(0,1,100)  
#H_begin = lambda f: 1.5+np.sin(2*np.pi*f)
H_begin = lambda f: (abs(f-0.05) < 0.01) + (abs(f-0.5)<0.005)

# define the matching frequencies: 
#  0 at the beginning remains zero after the warping
#  1 at the beginning remains one after warping
#  frequency 0.5 should be warped to 0.25 
f_matches = [(0,0), (1,1), (0.5,0.25)]

# generate the according inputs and outputs of the 2d function 
# that should be interpolated
f_in = np.array([x[0] for x in f_matches]*2)
t_in = np.array([0]*len(f_matches)+[1]*len(f_matches))
f_out = np.array([x[0] for x in f_matches] + [x[1] for x in f_matches])

# just printout the supporting points we have defined
for n in range(len(f_in)):
    print ("X(f=%f,t=%f)=%f" % (f_in[n], t_in[n], f_out[n]))

# generate the interpolation     
X = scipy.interpolate.interp2d(f_out, t_in, f_in)

# visualize the interpolation
t_range = np.linspace(0,1,100)
img = np.zeros((len(f), len(t_range)))
for n, t0 in enumerate(t_range):
    img[:,n] = H_begin(X(f,t0))
plt.imshow(img)
plt.xlabel('interpolation time')
plt.ylabel('frequency');

Output:

X(f=0.000000,t=0.000000)=0.000000
X(f=1.000000,t=0.000000)=1.000000
X(f=0.500000,t=0.000000)=0.500000
X(f=0.000000,t=1.000000)=0.000000
X(f=1.000000,t=1.000000)=1.000000
X(f=0.500000,t=1.000000)=0.250000

It does the following:

Essentially, you have a 2D-interpolation problem: You have pairs of frequencies that should be warped from the first element in the pair to the second element, when the interpolation factor $t$ increases from 0 to 1. (I.e. for $t=0$ you get the initial spectrum, for $t=1$ you get the final spectrum). So, given a mapping between some supporting start-frequencies and end-frequencies, we can create a 2-dimensional interpolation. For an input frequency and interpolation time $t$, it returns to where this frequency would be be warped to.

What is not solved:

• How to find the supporting points f_match. You might look for peak-matching algorithms or come up with something that is adapted for your application
• I think the interpolation should be done in log-domain rather than in linear domain if this is about audio processing. However, the algorithm wont change, just take the log of every frequency before interpolation and convert it back to linear frequency after interpolation.
• my script assumes you have accurate information on the spectrum at all frequencies (H_begin is a function). However, in reality you only have access to discrete frequencies of the spectrum. So, if the interpolation requires you to get the value from a frequency between two bins, you'd need to interpolate the original spectrum to approximate this value.