Fourier -- Sage

def conj(v): return vector(map(lambda z:z.conjugate(),v)) # Definition du conjugé d'un vecteur complexe

def produit_scalaire(v,w): return v*conj(w) # Definition du produit scalaire de vecteurs complexes - copier tel quel...

C=ComplexField(128) # truc à ne pas changer

N=16 # Nombre d'échantillons

B=list(vector(list(C(e^(i*2*pi*n*k/N)) for k in (0..N-1))) for n in (0..N-1)) # Définition de la base sur laquelle on projette - les "harmoniques". Ne pas changer!

plot(vector([imaginary(B[2][k]) for k in (0..N-1)])) # Un des vecteurs de la base...

Signal1=vector([RR(e^(-k/4)) for k in (0..N-1)]) # Signal 1 - à changer

plot(Signal1)

SpectreSignal1=vector([produit_scalaire(Signal1, B[n]) for n in (0..N-1)])

AmplitudesSignal1=vector(list(abs(SpectreSignal1[k]) for k in (0..N-1)))

AmplitudesSignal1

\newcommand{\Bold}[1]{\mathbf{#1}}\left(4.43801011026256,2.39867637710659,1.38122992711709,0.976588948128759,0.774510595394556,0.661460867672238,0.596559620922810,0.562519786080795,0.551879879103842,0.562519786080795,0.596559620922810,0.661460867672238,0.774510595394556,0.976588948128759,1.38122992711709,2.39867637710659\right)

show(plot(AmplitudesSignal1),ymin=0) # Changement pour avoir l'axe des x au bon endroit

DephasageSignal1=vector([arg(SpectreSignal1[k]) for k in (0..N-1)])

show(plot(DephasageSignal1))

SignalReconstruit1=(1/N)*sum(SpectreSignal1[k]*B[k] for k in (0..N-1))

SignalReconstruit1

\newcommand{\Bold}[1]{\mathbf{#1}}\left(1.00000000000000 + 2.77555756156289 \times 10^{-17}i,0.778800783071405,0.606530659712634 + 6.93889390390723 \times 10^{-18}i,0.472366552741014 + 6.93889390390723 \times 10^{-18}i,0.367879441171442 + 1.38777878078145 \times 10^{-17}i,0.286504796860190 - 2.77555756156289 \times 10^{-17}i,0.223130160148430 - 2.77555756156289 \times 10^{-17}i,0.173773943450445,0.135335283236613 + 1.38777878078145 \times 10^{-17}i,0.105399224561864 + 1.38777878078145 \times 10^{-17}i,0.0820849986238987 + 6.93889390390723 \times 10^{-18}i,0.0639278612067075,0.0497870683678639,0.0387742078317219,0.0301973834223185,0.0235177458560090 + 2.77555756156289 \times 10^{-17}i\right)

plot(vector([real(SignalReconstruit1[k]) for k in (0..N-1)]))

filtre1=vector( [ 1 for k in (0..7)] + [ 0 for k in (8..N-1)])

filtre1

\newcommand{\Bold}[1]{\mathbf{#1}}\left(1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0\right)

SpectreFiltre1=vector(list( filtre1[k]*SpectreSignal1[k] for k in (0..N-1)))

AmplitudesSignal1Filtre=vector(list(abs(SpectreFiltre1[k]) for k in (0..N-1)))

plot(AmplitudesSignal1Filtre)

SignalFiltre1=(1/N)*sum(SpectreFiltre1[k]*B[k] for k in (0..N-1))

plot(vector(list(real(SignalFiltre1[k]) for k in (0..N-1))))

def produit_convolution(v,w): return vector([sum(v[k]*w[n-k] for k in (0..N-1)) for n in (0..N-1)])

delta=vector([1]+[0 for k in (0..N-2)])

delta

\newcommand{\Bold}[1]{\mathbf{#1}}\left(1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0\right)

produit_convolution(Signal1,delta)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(1.00000000000000,0.778800783071405,0.606530659712633,0.472366552741015,0.367879441171442,0.286504796860190,0.223130160148430,0.173773943450445,0.135335283236613,0.105399224561864,0.0820849986238988,0.0639278612067076,0.0497870683678640,0.0387742078317220,0.0301973834223185,0.0235177458560091\right)

plot(_)

Spectre1delta=vector([ produit_scalaire(delta,B[n]) for n in (0..N-1)])

plot(vector([abs(Spectre1delta[n]) for n in (0..N-1)]) )

FiltreInverse1=(1/N)*sum(filtre1[k]*B[k] for k in (0..N-1))

plot(vector( [abs(FiltreInverse1[k]) for k in (0..N-1)]))

SignalConvolution=produit_convolution(FiltreInverse1,Signal1)

SpectreConvolution1=vector([produit_scalaire(SignalConvolution,B[k]) for k in (0..N-1)])

plot(vector([abs(SpectreConvolution1[k]) for k in (0..N-1)]))

Fourier

151 days ago by zigong