In [1]:
import control
import numpy as np

import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (10, 5)

def power_spectrum(t, d, fmin=0, fmax=30, width=0.5, n=8):
    T = t[-1] - t[0]
    ak = np.fft.ifft(d);
    dt = T/len(d)
    k = np.arange(0, len(d)//n)
    freq = k/T
    plt.bar(freq, np.abs(ak[k]), width=width)
    plt.ylabel('|ak|')
    plt.xlabel('Hz')
    plt.grid()
    plt.title('disturbance power spectrum')
    ax = plt.gca().set_xlim(fmin, fmax)
    return ak, freq
In [2]:
def gen_noise(filename, amps, freqs, phases, plot=False):
    dt = 0.001
    tf = 5
    t = np.arange(0, tf, dt)
    d = np.random.randn(len(t))*dt**(0.5)
    for amp, freq, phase in zip(amps, freqs, phases):
        d += amp*np.cos(2*np.pi*freq*t - phase)
    np.savetxt(filename, np.vstack([t, d]).T, delimiter=',', header='t, d',  fmt='%g')
    
    if plot:
        plt.figure()
        plt.plot(t, d)
        plt.grid()
        plt.xlabel('t, sec')
        plt.ylabel('y')
        plt.title(filename)

        plt.figure()
        power_spectrum(t, d)
        plt.title(filename + ' power_spectrum')