from numpy import * from matplotlib.pyplot import * # CONVOLUTION def convolute(A, B, T): N = max(len(A), len(B)) Ahat = fft.rfft(A,N) Bhat = fft.rfft(B,N) n = arange(N/2+1) prod = ((-1.0)**n) * Ahat * Bhat convoluted = T/N * fft.irfft(prod, N) return convoluted def crosscorr(A, B): ftA = fft.rfft(A) ftB = fft.rfft(B) return 1./len(A)*fft.irfft(ftA.conj()*ftB) N = 1000 T = 10.0 x = linspace(0, T, N) f = zeros(N) f[N*4.5/T:N*5.5/T] = 1.0 g = f g = zeros(N) g[N*1/T:N*3/T] = 0.5 g[N*6/T:N*7/T] = -1.0 figure(figsize=(8,4)) subplot(121) plot(x, f, "r-", linewidth=1.0) xlabel("\$x\$") ylabel("\$f(x)\$") axis([0,T,-1.2,1.2]) subplot(122) plot(x, g, "r-", linewidth=1.0) xlabel("\$x\$") ylabel("\$g(x)\$") axis([0,T,-1.2,1.2]) tight_layout() savefig("task_51.pdf") convolution = convolute(f, g, T) cc = T*crosscorr(f, g) figure(figsize=(8,4)) subplot(121) plot(x, convolution, "r-", linewidth=1.0) xlabel("\$x\$") ylabel("(\$f \star g)(x)\$") axis([0,T,-1.2,1.2]) subplot(122) plot(x, cc, "r-", linewidth=1.0) xlabel("\$x\$") ylabel("\$C(f,g)(x)\$") axis([0,T,-1.2,1.2]) tight_layout() show()