def generateSweepIR(self):
		Lsig=len(self.signal)
		Lavg=len(self.average[0,:])
		NFFTsig=nextpow2(Lsig)
		NFFTavg=nextpow2(Lavg)
		if NFFTsig != NFFTavg:
			raise Exception('NFFTsig != NFFTavg')
		else:
			NFFT=int(NFFTsig)
		sigfft=fft(self.signal,n=NFFT)
		print(self.average.shape)
		N=self.average.shape[0]
		avgfft=fftn(self.average,s=[NFFT])
		imp=np.array(np.zeros(avgfft.shape))
		#maximum=0
		for i in range(0,N):
			#imp[i,:]=np.real(fftshift(ifft(avgfft[i,:]/sigfft)))
			imp[i,:]=np.real(ifft(avgfft[i,:]/sigfft))
			imp[i,:]=fftshift(avgfft[:,i]/sigfft)
			imp[i,:find_nearest(f,-20000)]=0
			imp[i,find_nearest(f,-20):find_nearest(f,20)]=0
			imp[i,find_nearest(f,20000):]=0
			imp[i,:]=np.real(ifft(ifftshift(impfft)))
		#	if max(imp[i,:])>maximum:
		#		maximum=max(imp[i,:])
		#imp=imp/maximum
		#imp=imp/np.max(imp)
		toSave=np.array(imp.transpose(), dtype=np.float32)
		impfile=self.getImpFilename()#self.filepath.get()+os.sep+self.prefix.get()+'_IR_'+self.counter.get()+'.wav'
		scipy.io.wavfile.write(impfile,int(self.fs),toSave)
		raw=[]
		imp=[]

    def f(A_t, A_w, dz=delta_z):
        
        f.w_exp = exp(-1j * delta_z/2 * w_op)
    
        if f.delta_z != dz:
            f.initF()
            f.delta_z = dz
        
        ## Dispersion (I pass)
        f.A_w = A_w * f.w_exp
        f.A_t = ifft(f.A_w) / dt        
        
        ## Constant potential term
        f.A_t = f.A_t * f.t_exp

        ## Nonlinear operator as intensity dependency
        if nlin != 0:
            f.A_t *= exp(-1j * delta_z * nlin * absolute(f.A_t)**2)
            
        ## Additional nonlinear terms as a function t_nl_op(A(t),dt,z)
        if t_nl_op != None:
            f.A_t *= exp(-1j * delta_z * t_nl_op(f.A_t, dt, z0+delta_z/2) )

        ## Apodization
        if apod:
            f.A_t *= apod_array

        f.A_w = fft(f.A_t) * dt
        
        ## Dispersion (II pass)
        f.A_w *= f.w_exp
        f.A_t = ifft(f.A_w) / dt
        
        return f.A_t[:], f.A_w[:]

    def InvPotentialVorticity(self,field, length=None):

        if length is None:
            length = 2*pi;

        N = shape(field)[0];

        k = array(range(N),dtype=complex128);
        k = concatenate((range(0,N/2),range(-N/2,0)));
        k *= (2*pi)/length;

        [KX, KY] = meshgrid(k,k);

        """ We are trying to solve d_yy(eta) + d_xx(eta) - eta = p
        Therefore, in Fourier space, it will become
        (-(kx^2 + ky^2) - 1)etaHat = pHat
        """
        delsq = -(KX*KX + KY*KY) - 1 ;
        delsq[0,0] = 1;

        tmp = fft(field,axis=0);
        tmp = fft(tmp,axis=1);

        tmp = tmp/delsq;

        tmp = ifft(tmp,axis=1);
        tmp = ifft(tmp,axis=0);

        return tmp.real;

  def test_whitenize_with_decorator_size_4_without_fft_window(self, strategy):
    @strategy(size=4, hop=4)
    def whitenize(blk):
      return cexp(phase(blk) * 1j)

    sig = Stream(0, 3, 4, 0) # fft([0, 3, 4, 0]) is [7, -4.-3.j, 1, -4.+3.j]
                             # fft([4, 0, 0, 3]) is [7,  4.+3.j, 1,  4.-3.j]
    data4003 = ifft([1,  (4+3j)/5, 1,  (4-3j)/5]) # From block [4, 0, 0, 3]
    data0340 = ifft([1, (-4-3j)/5, 1, (-4+3j)/5]) # From block [0, 3, 4, 0]

    result = whitenize(sig) # No overlap-add window (default behavior)
    assert isinstance(result, Stream)

    expected = Stream(*data0340)
    assert almost_eq(result.take(64), expected.take(64))

    # Using a "triangle" as the overlap-add window
    wnd = [0.5, 1, 1, 0.5] # Normalized triangle
    new_result = whitenize(sig, ola_wnd=[1, 2, 2, 1])
    assert isinstance(result, Stream)
    new_expected = Stream(*data0340) * Stream(*wnd)
    assert almost_eq(new_result.take(64), new_expected.take(64))

    # With real overlap
    wnd_hop2 = [1/3, 2/3, 2/3, 1/3] # Normalized triangle for the new hop
    overlap_result = whitenize(sig, hop=2, ola_wnd=[1, 2, 2, 1])
    assert isinstance(result, Stream)
    overlap_expected = Stream(*data0340) * Stream(*wnd_hop2) \
                     + chain([0, 0], Stream(*data4003) * Stream(*wnd_hop2))
    assert almost_eq(overlap_result.take(64), overlap_expected.take(64))

def STRW(x, om, OM, up):
    M = len(x)
    N = len(x[0])
   
    # 2d convolution:
    y = np.ndarray(shape=(M, N), dtype='complex')
    h = get_wavelet(om, OM, up, M, N)
    
    # temporal convolution:
    for m in range(M):
        X = fft(x[m])
        H = fft(h[m])
        Z = X*H

        z = ifft(Z)
        y[m] = z

    # spectral convolution
    x = y.transpose()
    h = h.transpose()
    y = y.transpose()
    for n in range(N):
        X = fft(x[n])
        H = fft(h[n])
        Z = X*H
        
        z = ifft(Z)
        y[n] = z

    y = y.transpose()

    return y

def fwgn_model(fm,fs,N):
    N = int(N)
    Nfft = 2**max(3,nextpow2(2*fm/fs*N))
    Nifft = math.ceil(Nfft*fs/(2*fm))

    doppler_coeff = np.sqrt(doppler_filter(fm,Nfft))

    CGI, CGQ = fft(randn(Nfft)), fft(randn(Nfft))
    f_CGI = CGI * doppler_coeff
    f_CGQ = CGQ * doppler_coeff

    del CGI, CGQ, doppler_coeff
    gc.collect()

    tzeros = np.zeros(abs(Nifft-Nfft))
    filtered_CGI = np.hstack((f_CGI[:Nfft//2], tzeros, f_CGI[Nfft//2:]))
    filtered_CGQ = np.hstack((f_CGQ[:Nfft//2], tzeros, f_CGQ[Nfft//2:]))

    del tzeros, f_CGI, f_CGQ
    gc.collect()

    hI, hQ = ifft(filtered_CGI), ifft(filtered_CGQ)

    del filtered_CGI, filtered_CGQ
    gc.collect()

    rayEnvelope = np.abs(np.abs(hI) + 1j * hQ)
    rayRMS = math.sqrt(np.mean(rayEnvelope[:N]**2))

    # h_{I}(t) - jh_{Q}(t)
    # Here we have the phase shift of pi/2 when multiplying the imaginary
    # portion by -1j
    h = (np.real(hI[:N]) - 1j * np.real(hQ[:N]))/rayRMS

    return h

def autocorr_fft(signal, axis = -1):
    """Return full autocorrelation along specified axis. Use fft
    for computation."""
    if N.ndim(signal) == 0:
        return signal
    elif signal.ndim == 1:
        n       = signal.shape[0]
        nfft    = int(2 ** nextpow2(2 * n - 1))
        lag     = n - 1
        a       = fft(signal, n = nfft, axis = -1)
        au      = ifft(a * N.conj(a), n = nfft, axis = -1)
        return N.require(N.concatenate((au[-lag:], au[:lag+1])), dtype = signal.dtype)
    elif signal.ndim == 2:
        n       = signal.shape[axis]
        lag     = n - 1
        nfft    = int(2 ** nextpow2(2 * n - 1))
        a       = fft(signal, n = nfft, axis = axis)
        au      = ifft(a * N.conj(a), n = nfft, axis = axis)
        if axis == 0:
            return N.require(N.concatenate( (au[-lag:], au[:lag+1]), axis = axis), \
                    dtype = signal.dtype)
        else:
            return N.require(N.concatenate( (au[:, -lag:], au[:, :lag+1]), 
                        axis = axis), dtype = signal.dtype)
    else:
        raise RuntimeError("rank >2 not supported yet")

    def magabs_cs():
        #b.line_plot("magabs_cs", a.frequency, a.MagAbs[:, 0])
        #b.line_plot("magabs_cs", a.frequency, a.MagAbs[:, 257])

        if 0:
            myifft=fft.ifft(a.Magcom[:,500])
            myifft[50:-50]=0.0
            #myifft[:20]=0.0
            #myifft[-20:]=0.0
            bg=fft.fft(myifft)
            filt=[]
            for n in range(len(a.yoko)):
                myifft=fft.ifft(a.Magcom[:,n])
                #b.line_plot("ifft", absolute(myifft))
                myifft[50:-50]=0.0
                #myifft[:20]=0.0
                #myifft[-20:]=0.0
                filt.append(absolute(fft.fft(myifft)-bg))
            b.colormesh("filt", a.frequency, a.yoko, filt)
        if 1:
            myifft=fft.ifft(mag[:,500])
            myifft[40:-40]=0.0
            myifft[:20]=0.0
            myifft[-20:]=0.0
            bg=fft.fft(myifft)
            filt=[]
            for n in range(len(yok)):
                myifft=fft.ifft(mag[:,n])
                #b.line_plot("ifft", absolute(myifft))
                myifft[50:-50]=0.0
                #myifft[:20]=0.0
                #myifft[-20:]=0.0
                filt.append(absolute(fft.fft(myifft)))
            b.colormesh("filt", frq, yok, filt)

def error_test(vel, epsilon):
    sigma = 7.0
    tmax = 4 * L / vel
    dt = ((L / 2) / kmax) / 18
    psi[:] = exact(x, vel, 0.0, sigma)  # exp(-1.0*(x*x/16.0)+complex(0,vel)*x)
    psil[:] = exact(xl, vel, 0.0, sigma)
    m0 = sqrt(abs(psi[npoints / 4 : -npoints / 4] ** 2).sum())

    prop = exp(complex(0, -dt / 2.0) * k * k)
    propl = exp(complex(0, -dt / 2.0) * kl * kl)
    fil = tdpsf(npoints, dx, epsilon)

    nsteps = int(tmax / dt) + 1
    err = 0.0
    for i in range(nsteps):
        psi[:] = dft.ifft(prop * dft.fft(psi))  # Calculate TDPSF approximation
        fil(psi)

        psil[:] = dft.ifft(propl * dft.fft(psil))  # Calculate large box solution

        rem = abs(
            psil[npointsl / 2 - npoints / 4 : npointsl / 2 + npoints / 4]
            - psi[npoints / 2 - npoints / 4 : npoints / 2 + npoints / 4]
        )

        local_err = sqrt((rem * rem).sum()) / m0
        if local_err > err:
            err = local_err
            # print str(v) +", " + str(i*dt) + ", " + str(err)
        # print "T="+str(i*dt) + ", M = " + str(sqrt(abs(psi*psi).sum()*dx)/m0)
    print "k = " + str(vel) + ", epsilon = " + str(epsilon) + ", err = " + str(err)
    return err

    def lms(self):
        global wf_ech, wf_ref
        M = int(2 ** np.ceil(np.log2(self.M)))
        u = util.enframe(np.append(wf_ref, np.zeros(M)), 2 * M, M)
        d = util.enframe(wf_ech, M, 0)
        uf = nf.fft(u, 2 * M, 1)
        W = nf.fft(self.w, 2 * M)

        wf_ech = np.array([], np.complex128)

        for ii in range(0, u.shape[0]):
            mu_bar = (110 * np.exp(-(ii) / self.estep) + 1) * self.mu
            yfi = W * uf[ii, :]  # 1x2M
            yi = nf.ifft(yfi, 2 * M)  # 1x2M
            yi = yi[-M:]  # 1xM

            ei = d[ii, :] - yi  # 1xM
            efi = nf.fft(np.append(np.zeros(M), ei), 2 * M)  # 1x2M

            uefi = nf.ifft(np.conj(uf[ii, :]) * efi, 2 * M)
            uefi[-M:] = 0
            uei = nf.fft(uefi, 2 * M)

            W = W + 2 * mu_bar * uei
            wf_ech = np.append(wf_ech, ei)

        self.w = nf.ifft(W, 2 * M)

def DiffusiveFilter(sc,T,rho,c,L,dt=0.001,BW=0.7):
    
    # from numpy import linspace,zeros,vstack,array,dot,conj,pi
    # from scipy.signal import gausspulse
    # from numpy.fft import rfft,ifft

    t=linspace(0,T,round(T/dt))
    X=rfft(gausspulse((t-0.25*T),sc,bw=BW))
    
    s=linspace(0.,1/(2*dt),len(X))
    s[0]=1e-6
    
    RT=TransmissionReflection(s,rho,c,L)

    Y=X*RT
    
    s[0]=0.
    
    Y[:,0] = Y[:,1]
    
    y=ifft(2*Y,axis=1,n=2*Y.shape[1]-1)
        
    x=ifft(2*X,n=2*len(X)-1)
    
    t=linspace(0,T,y.shape[1])
    
    return t,vstack((x,y)),s,vstack((X,Y)),RT

def FFT_Correlation(x,y):
    """
    FFT-based correlation, much faster than numpy autocorr.
    x and y are row-based vectors of arbitrary lengths.
    This is a vectorized implementation of O(N*log(N)) flops.
    """

    lengthx = x.shape[0]
    lengthy = y.shape[0]

    x = np.reshape(x,(1,lengthx))
    y = np.reshape(y,(1,lengthy))

    length = np.array([lengthx, lengthy]).min()
    
    x = x[:length]
    y = y[:length]
    
    fftx = fft(x, 2 * length - 1, axis=1) #pad with zeros
    ffty = fft(y, 2 * length - 1, axis=1)

    corr_xy = fft.ifft(fftx * np.conjugate(ffty), axis=1)
    corr_xy = np.real(fft.fftshift(corr_xy, axes=1)) #should be no imaginary part

    corr_yx = fft.ifft(ffty * np.conjugate(fftx), axis=1)
    corr_yx = np.real(fft.fftshift(corr_yx, axes=1))

    corr = 0.5 * (corr_xy[:,length:] + corr_yx[:,length:]) / range(1,length)[::-1]
    return np.reshape(corr,corr.shape[1])

def _direct_fft_conv(log_pmf1, pmf1, fft1, log_pmf2, true_conv_len, fft_conv_len,
                     alpha, delta):
    if log_pmf2 is None:
        norms = np.linalg.norm(pmf1)**2
        fft_conv = fft1**2
    else:
        pmf2 = np.exp(log_pmf2)
        fft2 = fft.fft(pmf2, n = fft_conv_len)
        norms = np.linalg.norm(pmf1) * np.linalg.norm(pmf2)
        fft_conv = fft1 * fft2

    raw_conv = np.abs(fft.ifft(fft_conv)[:true_conv_len])
    error_level = utils.error_threshold_factor(fft_conv_len) * norms
    threshold = error_level * (alpha + 1)
    places_of_interest = raw_conv <= threshold

    if delta > error_level:
        ignore_level = delta - error_level

        conv_to_ignore = raw_conv < ignore_level
        raw_conv[conv_to_ignore] = 0.0
        places_of_interest &= ~conv_to_ignore

    log_conv = np.log(raw_conv)

    # the convolution is totally right, already: no need to consult
    # the support
    if not places_of_interest.any():
        return log_conv, []

    Q = 2 * len(log_pmf1) - 1 if log_pmf2 is None else len(log_pmf1) + len(log_pmf2) - 1

    support1 = log_pmf1 > NEG_INF
    if log_pmf2 is None:
        support2 = np.array([])
    else:
        support2 = log_pmf2 > NEG_INF
    if Q <= 2**42 and (not support1.all() or not support2.all()):
        # quickly compute the support; this is accurate given the
        # bound (if the vectors have no zeros, then the convolution
        # has no zeros, so we can skip this)
        fft_support1 = fft.fft(support1, n = fft_conv_len)
        if log_pmf2 is None:
            fft_support = fft_support1**2
        else:
            support2 = log_pmf2 > NEG_INF
            fft_support2 = fft.fft(support2, n = fft_conv_len)
            fft_support = fft_support1 * fft_support2

        raw_support = np.abs(fft.ifft(fft_support)[:true_conv_len])
        threshold = utils.error_threshold_factor(fft_conv_len) * 2 * Q
        zeros = raw_support <= threshold
        log_conv[zeros] = NEG_INF
        bad_places = np.where(~zeros & places_of_interest)[0]
    else:
        # can't/don't need to compute the support accurately.
        bad_places = np.where(places_of_interest)[0]

    return log_conv, bad_places

    def backward (self, f):

        from numpy import exp
        from numpy.fft import ifft

        f = ifft (fk, axis=1) 
        f *= exp (-1j*St*self.phase)
        return ifft (f, axis=0)

 def ifft_plot(self):
     p, pf=line(absolute(fft.ifft(self.Magcom[:,self.on_res_ind])), plotter="ifft_{}".format(self.name),
            plot_name="onres_{}".format(self.on_res_ind),label="i {}".format(self.on_res_ind), color="red")
     line(absolute(fft.ifft(self.Magcom[:,self.start_ind])), plotter=p, linewidth=1.0,
          plot_name="strt {}".format(self.start_ind), label="i {}".format(self.start_ind))
     line(absolute(fft.ifft(self.Magcom[:,self.stop_ind])), plotter=p, linewidth=1.0,
          plot_name="stop {}".format(self.stop_ind), label="i {}".format(self.stop_ind))
     return p

def ifft_plot(self):
    Magcom=(self.Magcom.transpose()-self.Magcom[:, 0]).transpose()
    p, pf=line(absolute(fft.ifft(Magcom[:,self.on_res_ind])), plotter="ifft_{}".format(self.name),
           plot_name="onres_{}".format(self.on_res_ind),label="i {}".format(self.on_res_ind))
    line(absolute(fft.ifft(Magcom[:,self.start_ind])), plotter=p,
         plot_name="strt {}".format(self.start_ind), label="i {}".format(self.start_ind))
    line(absolute(fft.ifft(Magcom[:,self.stop_ind])), plotter=p,
         plot_name="stop {}".format(self.stop_ind), label="i {}".format(self.stop_ind))

def OFDM_tx(IQ_data, Nf, N, Np=0, cp=False, Ncp=0):
    """
    Parameters
    ----------
    IQ_data : +/-1, +/-3, etc complex QAM symbol sample inputs
    Nf : number of filled carriers, must be even and Nf < N
    N : total number of carriers; generally a power 2, e.g., 64, 1024, etc
    Np : Period of pilot code blocks; 0 <=> no pilots
    cp : False/True <=> bypass cp insertion entirely if False
    Ncp : the length of the cyclic prefix

    Returns
    -------
    x_out : complex baseband OFDM waveform output after P/S and CP insertion

    Examples
    --------
    >>> import matplotlib.pyplot as plt
    >>> from sk_dsp_comm import digitalcom as dc
    >>> x1,b1,IQ_data1 = dc.QAM_bb(50000,1,'16qam')
    >>> x_out = dc.OFDM_tx(IQ_data1,32,64)
    >>> plt.psd(x_out,2**10,1);
    >>> plt.xlabel(r'Normalized Frequency ($\omega/(2\pi)=f/f_s$)')
    >>> plt.ylim([-40,0])
    >>> plt.xlim([-.5,.5])
    >>> plt.show()
    """
    N_symb = len(IQ_data)
    N_OFDM = N_symb // Nf
    IQ_data = IQ_data[:N_OFDM * Nf]
    IQ_s2p = np.reshape(IQ_data, (N_OFDM, Nf))  # carrier symbols by column
    print(IQ_s2p.shape)
    if Np > 0:
        IQ_s2p = mux_pilot_blocks(IQ_s2p, Np)
        N_OFDM = IQ_s2p.shape[0]
        print(IQ_s2p.shape)
    if cp:
        x_out = np.zeros(N_OFDM * (N + Ncp), dtype=np.complex128)
    else:
        x_out = np.zeros(N_OFDM * N, dtype=np.complex128)
    for k in range(N_OFDM):
        buff = np.zeros(N, dtype=np.complex128)
        for n in range(-Nf // 2, Nf // 2 + 1):
            if n == 0:  # Modulate carrier f = 0
                buff[0] = 0  # This can be a pilot carrier
            elif n > 0:  # Modulate carriers f = 1:Nf/2
                buff[n] = IQ_s2p[k, n - 1]
            else:  # Modulate carriers f = -Nf/2:-1
                buff[N + n] = IQ_s2p[k, Nf + n]
        if cp:
            # With cyclic prefix
            x_out_buff = fft.ifft(buff)
            x_out[k * (N + Ncp):(k + 1) * (N + Ncp)] = np.concatenate((x_out_buff[N - Ncp:],
                                                                       x_out_buff))
        else:
            # No cyclic prefix included
            x_out[k * N:(k + 1) * N] = fft.ifft(buff)
    return x_out

def D(uhat,kx,N,filt):
    """
    Function that computes the non linear term for the Burgers
    equation given the spectra of the variable. The dealiasing cuts
    the half of the spectra.
    """
    u = fft.ifft(uhat)
    dudx = fft.ifft(1.0j * kx * uhat)
    return fft.fft(u*dudx)*filt

def ComputeResponse(sc,T,rho,c,alpha,d):
    
    # from numpy import tan,pi,linspace,exp,zeros,hstack,vstack,array
    # from scipy.signal import gausspulse
    # from numpy.fft import rfft,ifft

    dt=1/(10*sc)
    t=linspace(0,T,round(T/dt))
    X=rfft(gausspulse((t-0.25*T),sc))
    
    Z=[]

    for i in range(len(rho)):
        Z.append(rho[i]*c[i])


    w=2*pi*linspace(0,1/(2*dt),len(X))

    R01=(Z[1]-Z[0])/(Z[1]+Z[0])
    T01=R01+1
    T10=1-R01
    R12=(Z[2]-Z[1])/(Z[1]+Z[2])
    T12=R12+1
    T21=1-R12

    Z234=Z[3]*(Z[4]-1j*Z[3]*tan(w*d[3]/c[3]))/(Z[3]-1j*Z[4]*tan(w*d[3]/c[3])) 
    R234=(Z234-Z[2])/(Z234+Z[2])

    T234=R234+1

    R45=(Z[5]-Z[4])/(Z[5]+Z[4])
    T45=R45+1
    

    Z432=Z[3]*(Z[2]-1j*Z[3]*tan(w*d[3]/c[3]))/(Z[3]-1j*Z[2]*tan(w*d[3]/c[3])) 
    R432=(Z432-Z[4])/(Z432+Z[4])

    T432=R432+1

    Y0=exp(-2*d[0]*alpha[0])*exp(-1j*w*2*d[0]/c[0])*R01*X
    Y1=exp(-2*d[1]*alpha[1])*exp(-1j*w*2*d[1]/c[1])*T01*T10*R12*Y0/R01

    Y2=exp(-2*d[2]*alpha[2])*exp(-1j*w*2*d[2]/c[2])*T12*T21*R234*Y1/R12
  
    Y3=exp(-2*d[4]*alpha[4])*exp(-1j*w*2*d[4]/c[4])*T234*T432*R45*Y2/R234
    
    Y4=exp(-2*d[4]*alpha[4])*exp(-1j*w*2*d[4]/c[4])*R432*R45*Y3

    Y5=exp(-2*d[4]*alpha[4])*exp(-1j*w*2*d[4]/c[4])*R432*R45*Y4

    y=2*ifft(vstack((Y0,Y1,Y2,Y3,Y4,Y5)),n=2*len(X)-1).transpose()

    x=2*ifft(X,n=2*len(X)-1)

    t=linspace(0,T,len(x))

    return t,x,y

def compute_solution():
    # frequency & start time
    f0 = f0_dominant_frequency
    t0 = t0_start_time

    # Size of the fourier transform
    N = 32*32*8
    T = T_length  # Length of the time signal
    t = linspace(0.,T,N)
    dt = t[1]-t[0]

    # source time function
    Rick = RickerF(t,t0,f0)
    # in frequency-domain
    Sf = f.fft(Rick)
    freq = f.fftfreq(N,d=dt)

    # plots figure
    fig = figure()
    ax = plot(freq[0:N/2],Sf[0:N/2],freq[0:N/2],abs(Sf[0:N/2]))

    #
    c0 = 1500
    k0 = 2*pi*f0/c0
    a0 = 0.000002

    print 'Q =',1/(c0*a0)

    # Position of the receiver
    PosX = receiver_X

    # Modification of the spectrum to include linear attenuation
    attf = k0*freq/f0 - 4 * a0 * freq * log((freq+0.0000001)/f0)
    Sf_att = Sf * exp(-1j * attf * PosX) * exp(-a0 * 2 * pi * freq * PosX)

    # Spectrum of the non attenuated signal
    Sf_0 = Sf * exp(-1j * 2 * pi * freq/c0 * PosX)

    #
    Sf_att_nan = nan_to_num(Sf_att)

    # plots figure
    plot(freq[0:N/2],Sf_att_nan[0:N/2],freq[0:N/2],abs(Sf_att_nan[0:N/2]))
    xlim(xmax=100)
    title('Spectra with and without attenuation')

    # signal in time-domain
    Sig_attn = f.ifft(Sf_att_nan)

    # plots figure
    figure(11)
    plot(t,Sig_attn)
    plot(t,f.ifft(Sf_0))
    grid()
    title('Snapshots with and without attenuation')
    show()