Untitled
raw download clone
TEXT
views 38
,
size 1830 b
import numpy as np
from math import *
import matplotlib.pyplot as plt

#Дискретное преобразование Фурье

def signal(x):
    return 1.0+sin(2.0*pi*x)+2.0*cos(4.0*pi*x)+0.5*cos(6.0*pi*x)

def dft(data):
    n=len(data)
    spectr=[0+0j]*n
    coef=2.0/n
    arg=pi*coef
    
    for j in range(n):
        spectr[j]=0.0+0.0j
        for i in range(n):
            spectr[j]+=data[i]*(cos(arg*i*j)-sin(arg*i*j)*1j)
    return np.array(spectr,dtype=complex)

#Обратное дискретное преобразование Фурье

def idft(data):
    n=len(data)
    signal=[0+0j]*n
    coef=2.0/n
    arg=pi*coef
    
    for j in range(n):
        signal[j]=0.0+0.0j
        for i in range(n):
            signal[j]+=data[i]*(cos(arg*i*(j+0))+sin(arg*i*(j+0))*1j)
    return np.array(signal,dtype=complex)/n

f=float(input('Опорная частота сигнала '))
T=float(input('Временной интервал '))
n=int(input('Число временных отчетов '))
t=np.linspace(0,T,n)
u=[0]*n

for i in range(n):
    u[i]=signal(f*t[i])
    
print('Расчет ДПФ:')
spec1=dft(u)
print('Расчет БПФ:')
spec2=np.fft.fft(u)
freq=np.fft.fftfreq(n,T/n)
 
print('Расчет OДПФ:')
sig_dft=idft(spec1)
print('Расчет OБПФ:')
sig_fft=np.fft.ifft(spec2)
 
plt.plot(t,u)
plt.show()
plt.plot(freq[0:n//2],(np.hypot(spec1.real,spec1.imag)/n*2.0)[0:n//2],label='ДПФ')
plt.plot(freq[0:n//2],(np.hypot(spec2.real,spec2.imag)/n*2.0)[0:n//2],'.-',label='БПФ')
plt.legend(loc='best')
plt.show()
plt.plot(t, sig_dft.real, label='ОДПФ')
plt.plot(t, sig_fft.real, '--', label='ОБПФ')
plt.plot(t,u, '--', label='Оригинал')
plt.legend(loc='best')
plt.axis(xmax=5)
plt.show()
close fullscreen
Login or Register to edit or fork this paste. It's free.