#!/usr/bin/env python
# coding: utf-8

# ## 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`](https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.random.random.html) 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]:


get_ipython().run_line_magic('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ę…