O co chodzi z tymi pikami¶

Wygenerujemy prostą funkcję okresową i obejrzymy jej transformatę Fouriera. Później będziemy ją zaburzać Funkcja np.random.random generuje wartości z przedziału $[0,1]$. Żeby otrzymać wartości z zakresu $[-1,1]$ odejmujemy od każdej wartości 0,5 i mnożymy ją przez 2.

In [1]:

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy

freq = 1 #Hz - liczba okresów na sekundę
amplitude = 3
time_to_plot = 100 # czas (długość danych wejściowych)
sample_rate = 20 # liczba próbek na sekundę
num_samples = sample_rate * time_to_plot

noise_amplitude = 0 # Zaczynamy od zerowych zaburzeń
noise = (np.random.random(num_samples) - 0.5) * 2
t = np.linspace(0, time_to_plot, num_samples)
sinus = [amplitude * np.sin(freq * i * 2*np.pi) for i in t]
signal = sinus + noise_amplitude * noise

Nasz sygnał wygląda tak:

In [2]:

plt.plot(t, signal);

In [3]:

fft_output = np.fft.rfft(signal)
magnitude_only = [np.sqrt(i.real**2 + i.imag**2)/len(fft_output) for i in fft_output]
frequencies = [(i*1.0/num_samples)*sample_rate for i in range(num_samples//2+1)]

In [4]:

plt.plot(frequencies, magnitude_only, 'r');

In [5]:

plt.plot(magnitude_only);
#plt.show()

Teraz zwiększamy amplitudę zaburzeń

In [6]:

noise_amplitude = 6
signal = sinus + noise_amplitude * noise
plt.plot(t, signal);

Na wykresie widać włąściwie tylko szum

In [7]:

fft_output = np.fft.rfft(signal)
magnitude_only = [np.sqrt(i.real**2 + i.imag**2)/len(fft_output) for i in fft_output]
frequencies = [(i*1.0/num_samples)*sample_rate for i in range(num_samples//2+1)]

In [8]:

plt.plot(frequencies, magnitude_only);

Mimo, że amplituda zaburzeń (6) jest dwukrotnie wyższa od amplitudy sygnału widać „pik” związany z częstością 1 Hz. Niestety, nie zawsze tak będzie

In [9]:

noise_amplitude = 14 # <---- cztery razy większa
signal = sinus + noise_amplitude * noise

fft_output = np.fft.rfft(signal)
magnitude_only = [np.sqrt(i.real**2 + i.imag**2)/len(fft_output) for i in fft_output]
frequencies = [(i*1.0/num_samples)*sample_rate for i in range(num_samples//2+1)]
plt.plot(frequencies, magnitude_only);

In [10]:

noise_amplitude = 15 # <---- pięć razy większa
signal = sinus + noise_amplitude * noise

fft_output = np.fft.rfft(signal)
magnitude_only = [np.sqrt(i.real**2 + i.imag**2)/len(fft_output) for i in fft_output]
frequencies = [(i*1.0/num_samples)*sample_rate for i in range(num_samples//2+1)]
plt.plot(frequencies, magnitude_only);

In [11]:

noise_amplitude = 18 # <---- sześć razy większa
signal = sinus + noise_amplitude * noise

fft_output = np.fft.rfft(signal)
magnitude_only = [np.sqrt(i.real**2 + i.imag**2)/len(fft_output) for i in fft_output]
frequencies = [(i*1.0/num_samples)*sample_rate for i in range(num_samples//2+1)]
plt.plot(frequencies, magnitude_only);

In [12]:

noise_amplitude = 24 # <---- osiem razy większa
signal = sinus + noise_amplitude * noise

fft_output = np.fft.rfft(signal)
magnitude_only = [np.sqrt(i.real**2 + i.imag**2)/len(fft_output) for i in fft_output]
frequencies = [(i*1.0/num_samples)*sample_rate for i in range(num_samples//2+1)]
plt.plot(frequencies, magnitude_only);

In [13]:

noise_amplitude = 30 # <---- osiem razy większa
signal = sinus + noise_amplitude * noise

fft_output = np.fft.rfft(signal)
magnitude_only = [np.sqrt(i.real**2 + i.imag**2)/len(fft_output) for i in fft_output]
frequencies = [(i*1.0/num_samples)*sample_rate for i in range(num_samples//2+1)]
plt.plot(frequencies, magnitude_only);

Zmniejszajac czas obserwacji sytuacja się pogarsza

In [14]:

freq = 1 #Hz - liczba okresów na sekundę
amplitude = 3
time_to_plot = 10 # czas (długość danych wejściowych) <--- !!!
sample_rate = 20 # liczba próbek na sekundę
num_samples = sample_rate * time_to_plot

noise_amplitude = 30
noise = (np.random.random(num_samples) - 0.5) * 2
t = np.linspace(0, time_to_plot, num_samples)
sinus = [amplitude * np.sin(freq * i * 2*np.pi) for i in t]
signal = sinus + noise_amplitude * noise

In [15]:

fft_output = np.fft.rfft(signal)
magnitude_only = [np.sqrt(i.real**2 + i.imag**2)/len(fft_output) for i in fft_output]
frequencies = [(i*1.0/num_samples)*sample_rate for i in range(num_samples//2+1)]
plt.plot(frequencies, magnitude_only, 'r');

Pik przy częstotliwości 1 widać, ale nie jest aż tak dominujący

In [16]:

freq = 1 #Hz - liczba okresów na sekundę
amplitude = 3
time_to_plot = 200 # czas (długość danych wejściowych) <--- !!!
sample_rate = 20 # liczba próbek na sekundę
num_samples = sample_rate * time_to_plot

noise_amplitude = 30
noise = (np.random.random(num_samples) - 0.5) * 2
t = np.linspace(0, time_to_plot, num_samples)
sinus = [amplitude * np.sin(freq * i * 2*np.pi) for i in t]
signal = sinus + noise_amplitude * noise
fft_output = np.fft.rfft(signal)
magnitude_only = [np.sqrt(i.real**2 + i.imag**2)/len(fft_output) for i in fft_output]
frequencies = [(i*1.0/num_samples)*sample_rate for i in range(num_samples//2+1)]
plt.plot(frequencies, magnitude_only, 'r');

Po zwiększeniu czasu obserwacji do 200 sekund — losowa struktura sygnału przebija się…