Fast Fourier Transform for Harmonic Analysis
Question: I'm analyzing the harmonics present in the wave motion as a function of where
along the string the pluck occurs. I hope to obtain a plot like those
exhibited on this page:
<https://softwaredevelopmentperestroika.wordpress.com/2013/12/10/fast-fourier-
transforms-with-python-the-noise-and-the-signal/>. To do this, I've written
code modelling an asymmetric triangle and implemented numpy's fft. The data
being output, however, is not what I expect, and it's peaked about a frequency
of 0 Hz. Here is my code and its output:

    
    
    from numpy.fft import fft as npfft, fftfreq as npfftfreq
    #triangular pulse
    def triangular_pulse(x, xmean, sigma):
        for i in x:
            if x[i]<=xmean:
                y[i] = x[i]*(sigma/xmean)
            else :
                y[i] = sigma-(x[i]-xmean)*(sigma/(200-xmean))
        return y
    
    N_masses = 200
    T  = 0.0669264714
    mu = .03937
    cSq = T/mu
    c  = np.sqrt(cSq)
    dx = 1.0
    
    dt = dx/c
    print dt
    
    #Initialize some arrays
    x0  = np.arange(N_masses)*dx
    y   = np.zeros(N_masses)
    vy  = np.zeros(N_masses)
    ay  = np.zeros(N_masses)
    
    #Set Initial conditions (pluck)
    # # half-pluck
    # y = 30*gaussian_pulse(x0,x0[N_masses/2],2)
    
    # quarter-pluck
    y = triangular_pulse(x0,x0[N_masses/4],1)
    
    rhat=npfft(y)
    freaq=npfftfreq(len(y),dt)
    plt.plot(freaq,np.abs(rhat)/len(rhat))
    plt.show()
    

Please let me know if you spot the source of my problem. Thanks!

# Update

Added y = triangular_pulse(x0,x0[N_masses/40],1) y-=np.mean(y) with the result
of a broader non-zero band; peak is still centered around "0", however.


Answer: Just subtract the mean of the signal before running the frequency analysis,
i.e. after calling `triangular_pulse`:

    
    
    y-=y.mean()
    

and you will obtain the peak at a non-zero frequency. This is because the
signal has a mean component that is not zero, that will show up as the
component at zero frequency.

EDIT: as a comment, you can rewrite the triangular pulse function using numpy
[where](http://docs.scipy.org/doc/numpy/reference/generated/numpy.where.html):

    
    
    def triangular_pulse2(x,xmean,sigma):
        return where(x<xmean,x*sigma/xmean,sigma-(x-xmean)*(sigma/(200-xmean)))
    

