Fourier Transformations for wave propagation: 

Introduction:  A Fast Fourier Transform converts wave form data in the time domain into the frequency domain. It achieves this by breaking down the original time-based waveform into a series of sinusoidal terms, each with a unique magnitude, frequency and phase. This process converts a waveform in the time domain into a series of sinusoidal functions which when added together reconstruct the original waveform. Plotting the amplitude of each sinusoidal term versus its frequency creates a power spectrum, which is the response of the original waveform in the frequency domain. 

 

When Fourier transforms are applicable, it means the “earth response” now is the same as the “earth response” later. Switching our point of view from time to space, the applicability of the Fourier transformation means that the “impulse response” here is the same as the “impulse response” there. An impulse is a column vector full of zeros with somewhere a one. An impulse response is a column from the matrix q = Bp The collection of impulse responses in q=Bp defines the convolution operation. 

 

Sample FFT application: 
import numpy as nm 

import scipy 

import scipy.fftpack 

import pylab 

 

def lowpass_cosine( y, tau, f_3db, width, padd_data=True): 

    # padd_data = True means we are going to symmetric copies of the data to the start and stop 

    # to reduce/eliminate the discontinuities at the start and stop of a dataset due to filtering 

    # 

    # False means we're going to have transients at the start and stop of the data 

 

    # kill the last data point if y has an odd length 

    if nm.mod(len(y),2): 

        y = y[0:-1] 

 

    # add the weird padd 

    # so, make a backwards copy of the data, then the data, then another backwards copy of the data 

    if padd_data: 

        y = nm.append( nm.append(nm.flipud(y),y) , nm.flipud(y) ) 

 

    # take the FFT 

    ffty=scipy.fftpack.fft(y) 

    ffty=scipy.fftpack.fftshift(ffty) 

 

    # make the companion frequency array 

    delta = 1.0/(len(y)*tau) 

    nyquist = 1.0/(2.0*tau) 

    freq = nm.arange(-nyquist,nyquist,delta) 

    # turn this into a positive frequency array 

    pos_freq = freq[(len(ffty)/2):] 

 

    # make the transfer function for the first half of the data 

    i_f_3db = min( nm.where(pos_freq >= f_3db)[0] ) 

    f_min = f_3db - (width/2.0) 

    i_f_min = min( nm.where(pos_freq >= f_min)[0] ) 

    f_max = f_3db + (width/2); 

    i_f_max = min( nm.where(pos_freq >= f_max)[0] ) 

 

    transfer_function = nm.zeros(len(y)/2) 

    transfer_function[0:i_f_min] = 1 

    transfer_function[i_f_min:i_f_max] = (1 + nm.sin(-nm.pi * ((freq[i_f_min:i_f_max] - freq[i_f_3db])/width)))/2.0 

    transfer_function[i_f_max:(len(freq)/2)] = 0 

 

    # symmetrize this to be [0 0 0 ... .8 .9 1 1 1 1 1 1 1 1 .9 .8 ... 0 0 0] to match the FFT 

    transfer_function = nm.append(nm.flipud(transfer_function),transfer_function) 

 

    # plot up the transfer function 

    # since "freq" is only the positive frequencies, select out 

    pylab.figure(1) 

    pylab.clf() 

    pylab.plot(freq,transfer_function) 

    pylab.xlabel('Frequency [Hz]') 

    pylab.ylabel('Filter Transfer Function') 

    pylab.xlim([-10.0,10.0]) 

    pylab.ylim([-0.05,1.05]) 

 

    # apply the filter, undo the fft shift, and invert the fft 

    filtered=nm.real(scipy.fftpack.ifft(scipy.fftpack.ifftshift(ffty*transfer_function))) 

 

    # remove the padd, if we applied it 

    if padd_data: 

        filtered = filtered[(len(y)/3):(2*(len(y)/3))] 

 

    # return the filtered data 

    return filtered 

 

 

# do an example of lowpass filtering 

# first make some fake data 

# a sine wave fluctuating once every pi seconds 

# samples 1000 times per second 

fakedata = nm.sin(nm.arange(0,11,0.001)) + nm.random.randn(len(nm.arange(0,11,0.001)))/4.0 

 

# run the filter 

# lowpass at 5 Hz, with a 1 Hz width of its roll-off 

filtered = lowpass_cosine(fakedata,0.001,5.0,1.0,padd_data=True) 

 

# plot the noisy data, with the filtered data on top 

pylab.figure(2) 

pylab.clf() 

pylab.plot(nm.arange(0,11,0.001),fakedata,label='Noisy Data') 

pylab.plot(nm.arange(0,11,0.001),filtered,label='Lowpass Filtered Data') 

pylab.xlabel('Time [s]') 

pylab.ylabel('Voltage') 

pylab.legend() 

 

pylab.ion() 

pylab.show()