-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathfft.py
202 lines (164 loc) · 6.32 KB
/
fft.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
# Written by serhatsoyer
import numpy as np
import matplotlib.pyplot as plt
def get_time(first_time, last_time, samp_freq):
"""
Generates a time vector in numpy
from first_time (s) to last_time (s) at samp_freq (Hz)
"""
return np.linspace(first_time, last_time, 1 + (samp_freq * (last_time - first_time)), endpoint=True)
def get_signal(frequency, amplitude, phase, offset, time):
"""
Generates a sinusoid at frequency (Hz)
amplitude is in signal units
phase is in radians
offset is in signal units
time vector is produced by get_time function
"""
return offset + amplitude * np.sin(2 * np.pi * frequency * time + phase)
def remove_offset(signal):
"""
Removes the mean of the signal generated by get_signal (signal)
"""
return signal - np.mean(signal)
def get_freq_range(signal, samp_freq):
"""
Gets the x-axis of the one-sided fft in Hz
signal is generated by get_signal
samp_freq is in Hz
"""
return samp_freq * np.fft.fftfreq(len(signal))[:int(len(signal) / 2)]
def take_fft(signal):
"""
Takes the fft of the signal and scales the output such that
the output unit is equal to the signal unit
The output fft is one-sided
signal is produced by get_signal
"""
fft = (abs(2 * np.fft.fft(signal) / len(signal)))[:int(len(signal) / 2)]
fft[0] /= 2
return fft
def get_total_power(time, signal, window_length=256):
"""
Estimates the power of windowed signal using np.std()
time is generated by get_time function
signal is generated by get_signal function
Returns both output time vector and the estimated power vector for each window
Time elements correspond to the middle points of the windows
"""
power_time = []
power_signal = []
half_window_length = round(window_length / 2)
for left_idx in range(min(len(time), len(signal)) - window_length):
power_time.append(time[left_idx + half_window_length])
power_signal.append(np.std(remove_offset(signal[left_idx:(left_idx + window_length)])))
return power_time, power_signal
def plot_signal_to_axis(axis, time, signal):
"""
Used by plot_signal, plot_signal_and_fft, and plot_signal_and_spectrogram
"""
axis.plot(time, signal)
axis.grid(True)
axis.set_xlabel('Time (s)')
axis.set_ylabel('Signal Unit')
def plot_fft_to_axis(axis, signal, samp_freq):
"""
Used by plot_fft and plot_signal_and_fft
"""
axis.stem(get_freq_range(signal, samp_freq), take_fft(signal))
axis.grid(True)
axis.set_xlabel('Frequency (Hz)')
axis.set_ylabel('Signal Unit')
def plot_spectrogram_to_axis(axis, time, signal, samp_freq):
"""
Used by plot_spectrogram, plot_signal_and_spectrogram, and plot_spectrogram_and_power
"""
axis.specgram(signal, Fs=samp_freq)
axis.set_xlabel('Time (s)')
axis.set_ylabel('Frequency (Hz)')
axis.set_xlim(time[0], time[-1])
def plot_total_power_to_axis(axis, time, signal):
"""
Used by plot_total_power and plot_spectrogram_and_power
"""
power_time, power_signal = get_total_power(time, signal)
axis.plot(power_time, power_signal)
axis.grid(True)
axis.set_xlabel('Time (s)')
axis.set_ylabel('Signal Unit')
def plot_signal(time, signal):
"""
Plots signal vs. time
time is generated by get_time
signal is generated by get_signal
"""
plt.style.use('Solarize_Light2')
default_size = plt.rcParams.get('figure.figsize')
_, axis = plt.subplots(1, 1, figsize=(default_size[0], default_size[1]))
axis.set_title('Signal in Time Domain')
plot_signal_to_axis(axis, time, signal)
def plot_fft(signal, samp_freq):
"""
Plots one-sided and nicely scaled fft of the signal at samp_freq Hz
signal is generated by get_signal
"""
plt.style.use('Solarize_Light2')
default_size = plt.rcParams.get('figure.figsize')
_, axis = plt.subplots(1, 1, figsize=(default_size[0], default_size[1]))
axis.set_title('Signal in Frequency Domain')
plot_fft_to_axis(axis, signal, samp_freq)
def plot_spectrogram(time, signal, samp_freq):
"""
Plots the spectrogram of the signal at samp_freq Hz
time is generated by get_time
signal is generated by get_signal
"""
plt.style.use('Solarize_Light2')
default_size = plt.rcParams.get('figure.figsize')
_, axis = plt.subplots(1, 1, figsize=(default_size[0], default_size[1]))
axis.set_title('Spectrogram')
plot_spectrogram_to_axis(axis, time, signal, samp_freq)
def plot_total_power(time, signal):
"""
Plots the estimated total power of the windowed signal
time is generated by get_time
signal is generated by get_signal
"""
plt.style.use('Solarize_Light2')
default_size = plt.rcParams.get('figure.figsize')
_, axis = plt.subplots(1, 1, figsize=(default_size[0], default_size[1]))
axis.set_title('Total Power')
plot_total_power_to_axis(axis, time, signal)
def plot_signal_and_fft(time, signal, samp_freq):
"""
Plots the signal and the fft of the signal in subplots
"""
plt.style.use('Solarize_Light2')
default_size = plt.rcParams.get('figure.figsize')
_, axes = plt.subplots(2, 1, figsize=(default_size[0], 1.5*default_size[1]))
axes[0].set_title('Signal in Time and Frequency Domains')
plot_signal_to_axis(axes[0], time, signal)
plot_fft_to_axis(axes[1], signal, samp_freq)
plt.show()
def plot_signal_and_spectrogram(time, signal, samp_freq):
"""
Plots the signal and the spectrogram of the signal in subplots
"""
plt.style.use('Solarize_Light2')
default_size = plt.rcParams.get('figure.figsize')
_, axes = plt.subplots(2, 1, figsize=(default_size[0], 1.5*default_size[1]))
axes[0].set_title('Signal and Spectrogram')
plot_signal_to_axis(axes[0], time, signal)
plot_spectrogram_to_axis(axes[1], time, signal, samp_freq)
plt.show()
def plot_spectrogram_and_power(time, signal, samp_freq):
"""
Plots spectrogram and the estimated total power of the signal in subplots
"""
plt.style.use('Solarize_Light2')
default_size = plt.rcParams.get('figure.figsize')
_, axes = plt.subplots(2, 1, figsize=(default_size[0], 1.5*default_size[1]))
axes[0].set_title('Spectrogram and Total Power')
plot_spectrogram_to_axis(axes[0], time, signal, samp_freq)
plot_total_power_to_axis(axes[1], time, signal)
plt.show()