## 1. Shannon's theorem in use.¶

In [ ]:
import numpy as np
import scipy.fftpack as sf
import matplotlib.pyplot as plt
import pylab as pl
from __future__ import division


Hereafter, we define a continuous signal that we sample and then we use Shannon's interpolation formula.

In [ ]:
def mySignal(t):
fq = 0.9
fq2 = 0.3
return np.sin(2*np.pi * fq*t) + 0.2*np.cos(2*np.pi * fq2*t)
x = np.linspace(0,2*np.pi,1000)
plt.plot(x,mySignal(x))
plt.show()


### Let us now observe the aliasing phenomenon by subsampling an image which presents high-frequency information.¶

In [ ]:
import imageio as imio
sub_defense = colored_image[::5,::5,:]
print(np.shape(colored_image))
#plt.gray()
plt.figure(figsize = (20,20))
plt.subplot(1,2,1)
plt.title("initial image")
plt.imshow(colored_image)
plt.subplot(1,2,2)
plt.title("sub-sampled image")
plt.imshow(sub_defense)

In [ ]:
import imageio as imio
defense = np.sum(colored_image*[ 0.21, 0.72 ,0.07],axis=-1)
sub_defense = defense[::2,::2]
print(np.shape(defense))
plt.gray()
plt.figure(figsize = (20,20))
plt.subplot(1,2,1)
plt.title("arche")
plt.imshow(defense)
plt.subplot(1,2,2)
plt.title("arche sub-sampled")
plt.imshow(sub_defense)


## A. On a 1D signal.¶

In [ ]:
# Choose a grid x and a signal y
n = 1000
m = int(n/2)
print(m)
x = np.linspace(-1,1,n)
y = np.zeros_like(x)
y[m:] = 1
y[0:m] = -1
plt.plot(x,y)
plt.title(" A step function ");

In [ ]:
# Use of FFT to compute the discrete Fourier transform
import scipy.fftpack as sf
spectre = sf.fft(y)
inverse_du_spectre = sf.ifft(spectre)


## Important remark:¶

The FFT output gives the DFT of the input for positive frequencies starting from $0$ on the first half of the vector and the negative frequencies on the second half starting always in increasing order. More precisely, we have that, if $y$ is the output, $[y(0),y(1),..,y(n/2),y(1-n/2),...,y(-1)]$ if $n$ is even, $[y(0),y(1),..,y((n-1)/2),y(-(n-1)/2),...,y(-1)]$ otherwise. In general, always use $fftshift$ to put low frequencies in the middle, which is the customary representation of the FFT.

In [ ]:
# illustration.
plt.figure(figsize = (10,6))
plt.subplot(1,2,1)
plt.plot(np.abs(spectre))
plt.title(" The spectrum ")
plt.subplot(1,2,2)
plt.plot(np.abs(sf.fftshift(spectre)))
plt.title(" The 0 frequency is in the middle. ")

In [ ]:
reduced = np.copy(spectre)
h = 40
reduced[h:n-h] = 0
reconstruct=sf.ifft(reduced)
plt.plot(np.real(reconstruct))


## Answer. $#####################################$¶

In [ ]:
plt.title("Decreasing of the spectrum log-modules")
fq = spectre[0:m]
plt.plot(np.log(np.abs(fq)));


## B. On an image.¶

In [ ]:
import scipy.signal as sig
from scipy import misc

In [ ]:
# load an image
image = misc.face(gray=True).astype(np.float32)
plt.figure(figsize = (10,10))
plt.subplot(1,2,1)
plt.gray()
plt.title(" Grey Level Image " + str(np.shape(image)))
plt.imshow(image)
# sub-sample a matrix
image_sub = image[::3,::4]
plt.subplot(1,2,2)
plt.gray()
plt.title(" Sub-sampled Image " + str(np.shape(image_sub)))
plt.imshow(image_sub)


### Hereafter, we represent the logarithm of the module of the image spectrum.¶

In [ ]:
spectre_im = sf.fft2(image_sub)
plt.figure(figsize = (10,10))
plt.subplot(1,2,1)
plt.gray()
plt.title(" Use of $fftshift$ in 2D ")
plt.imshow(sf.fftshift(image_sub))
plt.subplot(1,2,2)
plt.gray()
plt.title(" Image spectrum ")
plt.imshow(sf.fftshift(np.log(np.abs(spectre_im))))


### Ex: Explain what is the reason of appearance of the white lines in the frequency representation of the image. Propose and implement a solution to get rid of this artifact. We advise to use a mask on the intial image defined using $\sin$ or $\cos$ and then to compute the FFT.¶

In [ ]:
# Your code here.

In [ ]:
# Your code here.


## Hereafter, we construct an image from a random modification of the phase of the FFT of an image, which contains what is called texture information, that is repetitive patterns.¶

In [ ]:
import imageio as imio
texture_im = np.sum(bois*[ 0.21, 0.72 ,0.07],axis=-1)
spectre_im = sf.fft2(texture_im)
# construct the random phase.
plt.figure(figsize = (10,10))
temp = 0.3 * np.random.rand(*np.shape(spectre_im))
random_phase = np.cos(2*np.pi * temp) + np.sin(2*np.pi*temp)*1j
new_spectrum = spectre_im*random_phase
temp_2 = 30 * temp
random_phase_2 = np.cos(2*np.pi * temp_2) + np.sin(2*np.pi*temp_2)*1j
new_spectrum_2 = spectre_im*random_phase_2
plt.subplot(1,3,1)
plt.title(" Initial image ")
plt.imshow(texture_im)
plt.subplot(1,3,2)
plt.title(" with a random phase ")
plt.imshow(np.real(ifft2(new_spectrum)))
plt.subplot(1,3,3)
plt.title(" Larger random phase ")
plt.imshow(np.real(ifft2(new_spectrum_2)))


## 3. Convolution.¶

In [ ]:
### Construct a gaussian filter.


## Ex: Suppressing the aliasing effect.¶

### Start with the arch image at full resolution, convolve it with a gaussian filter and subsample it as done in the aliasing experiment. Vary the constant of the gaussian filter and try to suppress the aliasing.¶

In [ ]:
### your code here.


## 3. Linear and non-linear approximation.¶

In [ ]:
### We define two measures to quantify the quality of the signal with respect to the ground truth.
### It is a quantity that is computed between x,y two images. They are called peaked signal to noise ratio,
### and signal to noise ration (psnr and snr).
def psnr(x, y, vmax=-1):
d = np.mean((x - y) ** 2)
if d ==0:
return "Equal inputs"
if vmax < 0:
m1 = abs(x).max()
m2 = abs(y).max()
vmax = max(m1, m2)

return 10 * np.log10(vmax ** 2 / d)

def snr(x, y):
s =  np.linalg.norm(x - y)
if s == 0:
return "Equal inputs"
return 20 * np.log10(np.linalg.norm(x) /s)

In [ ]:
def LinearApproximation(x,n):
if 2*n>=np.max(np.shape(x)):
print "n out of bounds"
return x
spectre = sf.fftshift(sf.fft2(x))
filtre = np.zeros_like(x)
filtre[n:-n,n:-n] = 1
result = np.real(sf.ifft2(sf.fftshift(filtre*spectre)))
return result
approx1 = LinearApproximation(image_sub,60)
approx2 = LinearApproximation(image_sub,115)
plt.figure(figsize = (15,15))
plt.subplot(1,3,1)
plt.title(" snr " + str(round(psnr(image_sub,approx1),3)))
plt.imshow(approx1)
plt.subplot(1,3,2)
plt.title(" snr " + str(round(psnr(image_sub,approx2),3)))
plt.imshow(approx2)
plt.subplot(1,3,3)
plt.title(" Initial image ")
plt.imshow(image_sub);


### Compare the linear and the non-linear approximation.¶

In [ ]:
### Your code here.


## 4. STFT (Short Time Fourier Transform) and source separation.¶

### In this part, we use STFT which is a collection of Fourier transform of a 1D signal on time subintervals. We use it in order to experiment source separation which comes from the idea that the different signals may present very different behaviour in the frequency domain.¶

In [ ]:
# utils to load the sounds.
import numpy as np
import wave as wv

x_raw = wv.open(file)
n = x_raw.getnframes()
x_raw.close()

if file[::-1][:8][::-1] == "bird.wav":
x = np.delete(x,list(range(6001)) + list(range(12500, 15001)) + list(range(22500, 24001)) + list(range(32500,34001)))

if n0 !=0 and n0 < n:
x = x[:n0]

return x/np.max(x)


## Hereafter, we load the $3$ sounds and plot the second one.¶

In [ ]:
n = 1024*16
s = 3 #number of signals.
x = np.zeros([n,3])

x = x/np.tile(np.std(x,0),(n,1))
p = 2 #number of micros
import matplotlib.pyplot as plt
plt.figure(figsize = (10,10))
plt.xlim(0,n)
plt.plot(x[:,2])


### We simulate two micros which are implemented by linear combinations of the signals.¶

In [ ]:
theta = np.linspace(0, np.pi, s + 1)[:-1]
theta[0] = .2
M = np.vstack((np.cos(theta), np.sin(theta)))
## recovered signals
y = np.dot(x,np.transpose(M))


### We use the STFT function from the python package signal and plot it.¶

In [ ]:
import scipy.signal as sig
f,t,w = sig.stft(x[:,0])
a,b,z = sig.stft(x[:,1])
plt.figure(figsize = (10,10))
plt.subplot(1,2,1)
plt.title("stft of signal 1")
plt.imshow(np.log(np.abs(w)))
plt.subplot(1,2,2)
plt.title("stft of signal 2")
plt.imshow(np.log(np.abs(z)))


### We numerically check that the STFT and ISTFT are indeed inverse from each others.¶

In [ ]:
micro1 = y[:,0]
micro2 = y[:,1]
stft = lambda im : sig.stft(im,noverlap = 64,nperseg = 128)
istft = lambda im : sig.istft(im,noverlap = 64,nperseg = 128)
f,t,w1 = stft(micro1)
f,t,w2 = stft(micro2)
t,recov = istft(w1)
print(np.sum((micro1 - recov)**2))


### By selecting randomly a group of points in time, we plot in the plane the coordinates of the points, being the measured signals by the micros.¶

In [ ]:
nbre_selec = 600
from random import shuffle
print(np.shape(y)[0])
liste = range(np.shape(y)[0])
shuffle(liste)
plot(y[liste[0:nbre_selec],0],y[liste[0:nbre_selec],1],"o",markersize=1)


### Q: Do the same with the STFT signals, what do you observe ?¶

In [ ]:
nbre_selec = 3000
from random import shuffle
H = np.asarray([np.imag(w1.flatten()),np.imag(w2.flatten())]).transpose()
liste2 = range(np.shape(H)[0])
shuffle(liste2)
plt.figure(figsize = (10,10))
plt.plot(H[liste2[0:nbre_selec],0],H[liste2[0:nbre_selec],1],"o",markersize=2);


### Q. Comment the code below after running it.¶

In [ ]:
#import math
#Theta = np.zeros(np.shape(H)[0])
#for i in range(np.shape(H)[0]):
#    Theta[i] = math.atan2(H[i,1],H[i,0])%np.pi
#print(np.shape(Theta))

In [ ]:
#nbins = 400
#t = np.linspace(np.pi/200,np.pi,nbins)
#hist = np.histogram(Theta[:10000],t)
#h = hist[0]/np.sum(hist[0])
#t = t[:-1]

#plt.figure(figsize = (7,5))
#plt.bar(t, h, width = np.pi/nbins, color = "darkblue", edgecolor = "darkblue")
#plt.xlim(0,np.pi)
#plt.ylim(0,np.max(h))
#plt.show()

In [ ]:
## Affiche les valeurs d'angles les plus importantes.
#theta = []
#for xx in (h>0.01)*t:
#    if xx>0:
#        theta.append(xx)
#print(theta)


### Q. Comment the code below after running it ?¶

In [ ]:
projections = np.dot(Estimated,W)
C = np.abs(projections)
temp = np.max(C,0)
I = np.argmax(C,0)
threshold = .005
D = np.sqrt(np.sum(W**2, 0))
print(np.shape(D))
I = I*(D > threshold)

In [ ]:
masque = np.zeros_like(projections)
for i in range(np.shape(masque)[0]):
masque[i] = 1*(I==i)
source_stft = projections * masque
print(np.shape(source_stft))

In [ ]:
MaListe = []
for i in range(3):
f,t,w = stft(x[:,i])
MaListe.append(w)
print(np.shape(w))
print(i,snr(np.abs(w.flatten()),np.abs(source_stft[i,:])))

In [ ]:
print(np.shape(MaListe[1]))

In [ ]:
i = 1
plt.figure(figsize = (40,40))
plt.subplot(1,2,1)
plt.title("stft of source")
plt.imshow(np.log(np.abs(MaListe[i])))
plt.subplot(1,2,2)
plt.title("stft of recovered signal")
plt.imshow(np.log(np.abs(source_stft[i,:]) + 1e-10).reshape(np.shape(MaListe[i])))

In [ ]:
X = []
for i in range(3):
temp = source_stft[i,:].reshape(np.shape(MaListe[i]))
X.append(istft(temp)[1])
h = X[0]
print(np.shape(h))

In [ ]:
i = 0
l = 600
m = 8000
plt.figure(figsize = (10,10))
plt.subplot(1,2,1)
plt.title("recovered")
plt.xlim(l,m)
plt.plot(X[i])
plt.subplot(1,2,2)
plt.title("source")
plt.xlim(l,m)
plt.plot(x[:,i])
print(snr(X[i],x[:,i]))