Ревущий шумовой сигнал fft - PullRequest
Новая
       80

Ревущий шумовой сигнал fft

0 голосов
/ 29 мая 2020

Я хотел создать два случайных шумовых сигнала с дискретизацией 2,5 Гвыб / с в диапазоне частот 200 кГц - 20 МГц, время сигнала 5 мкс и вычислить его fft, но у меня проблема с fft. Спасибо за помощь, код 

import numpy as np
import matplotlib.pyplot as plot
from scipy import signal
from scipy import fft
import pandas as pd


   t = np.arange(0, 5e-6, 4e-10)

   s1 = 1e-8*np.random.normal(0, 1, 12500)
   s2 = 1e-8*np.random.normal(0, 1, 12500)
   sos1 = signal.butter(N=10, Wn=[200000, 20000000], btype='band', fs=2.5e9, output='sos')
   sos2 = signal.butter(N=10, Wn=[200000, 20000000], btype='band', fs=2.5e9, output='sos')

   fs1 = signal.sosfilt(sos1, s1)
   fs2 = signal.sosfilt(sos2, s2)

f1 = abs(fs1.fft())
f2 = abs(fs2.fft())
ax1 = plot.subplot(311)
plot.plot(t, fs1, t, fs2)
#ax1.set_xlim([0, 5e-6])
plot.xlabel('Time (s)')
plot.ylabel('Current (A)')
ax2 = plot.subplot(312)
plot.plot(f1, f2)
plot.xlabel('Frequency (Hz)')
plot.ylabel('Current (A)')

plot.show()

1 Ответ

1 голос
/ 29 мая 2020

Мне пришлось внести некоторые изменения в ваш код, чтобы запустить его. Основная проблема заключалась в замене fs1.fft() на fft.fft().

Другой проблемой является метод 'fft.fftshift ()', о котором вам следует знать. Вы можете вычислить частотный вектор вручную, но это несколько утомительно из-за порядка элементов в результирующем векторе fft. Результат fft имеет своеобразное расположение частот. Из документации scipy.fft.fft () :

Частотный член f = k / n находится в y [k]. При y [n / 2] мы достигаем частоты Найквиста и переходим к членам с отрицательной частотой. Итак, для 8-точечного преобразования частоты результата равны [0, 1, 2, 3, -4, -3, -2, -1]. Чтобы изменить порядок вывода fft так, чтобы компонент нулевой частоты был центрирован, например [-4, -3, -2, -1, 0, 1, 2, 3], используйте fftshift.

Итак, самый простой способ - использовать scipy.fft.fftfreq () , чтобы scipy сделал расчет за вас. Если вы хотите построить его естественным образом, вам следует вызвать scipy.fft.fftshift () , чтобы сдвинуть частоту нулевой Гц в центр массива.

Также, как вы используете реальные сигналы, из соображений эффективности вы можете рассмотреть возможность использования вещественной версии алгоритма fft scipy.fft.rfft () . Выходные данные не включают отрицательные частоты, поскольку для реальных входных массивов выход полного алгоритма всегда симметричен c.

См. Код ниже. 

import matplotlib
matplotlib.use('Qt5Agg')

import numpy as np
import matplotlib.pyplot as plot
from scipy import signal
from scipy import fft
import pandas as pd

sampling_freq_Hz = 2.5e9
sampling_period_s = 1 / sampling_freq_Hz
signal_duration_s = 5.0e-6
wanted_number_of_points = signal_duration_s / sampling_period_s
f_low_Hz = 200.0e3
f_high_Hz = 20.0e6
msg = f'''
Sampling frequency: {sampling_freq_Hz} Hz
Sampling period: {sampling_period_s} s
Signal duration: {signal_duration_s} s
Wanted number of points: {wanted_number_of_points}
Lower frequency limit {f_low_Hz}
Upper frequency limit: {f_high_Hz}
'''
print(msg)

# Time axis
time_s = np.arange(0, signal_duration_s, sampling_period_s)
real_number_of_points = time_s.size
print(f'Real number of points: {real_number_of_points}')

# Normal(0,sigma^2) distribuited random noise
sigma_2 = 1.0e-8
s1 = np.random.normal(0, sigma_2, real_number_of_points)
s2 = np.random.normal(0, sigma_2, real_number_of_points)

# Since both filters are equal, you only need one 
sos1 = signal.butter(N=10, Wn=[f_low_Hz, f_high_Hz], btype='band', fs=sampling_freq_Hz, output='sos')
#sos2 = signal.butter(N=10, Wn=[f_low_Hz, f_high_Hz], btype='band', fs=sampling_freq_Hz, output='sos')

# Do the actual filtering
filtered_signal_1 = signal.sosfilt(sos1, s1)
filtered_signal_2 = signal.sosfilt(sos1, s2)

# Absolute value
f_1 = abs(fft.fft(filtered_signal_1))
f_2 = abs(fft.fft(filtered_signal_2))
freqs_Hz = fft.fftfreq(time_s.size, sampling_period_s)

# Shift the FFT for understandable plotting
f_1_shift = fft.fftshift(f_1)
f_2_shift = fft.fftshift(f_2)
freqs_Hz_shift = fft.fftshift(freqs_Hz)

# Plot
ax1 = plot.subplot(311)
ax1.plot(time_s, filtered_signal_1, time_s, filtered_signal_2)
#ax1.set_xlim([0, 5e-6])
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Current (A)')

ax2 = plot.subplot(313)
ax2.plot(freqs_Hz_shift, f_1_shift, freqs_Hz_shift, f_2_shift)
ax2.set_xlabel('Frequency (Hz)')
ax2.set_ylabel('Current (A)')

plot.show()
...