Contents

TD2 - FFT

1. Shannon’s theorem in use.

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

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

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()
../../_images/TP2 - FFT - Exercices_4_0.png

Ex. Choose a sampling time \(T_e\), sample this signal and then show the reconstruction using linear interpolation.

fq = 0.9
fq2 = 0.3
Te = 1/(2*max(fq,fq2))*0.9
# fe = 1/Te
N = int((2*np.pi)/Te)
sampled_signal = []
x0 = []

for i in range(0, N-1):
    n = i * Te
    sampled_signal.append(mySignal(n))
    x0.append(n)

pl.plot(x0, sampled_signal)
pl.show()
../../_images/TP2 - FFT - Exercices_6_0.png

Ex. Implement Shannon’s reconstruction function and implement the associated reconstruction formula. Change the input signal and compare with linear interpolation. What do you observe?

\(sinc(t) = \frac{sin(\pi t)}{(\pi t)} \) avec \(sinc(0) = 1\)

def sinc(t):
    x_temp = (np.pi*t) / Te
    return 1 if t == 0 else np.sin(x_temp) / x_temp

def shannon_reconstruct(t):
    signal_sum = 0
    for i in range(0, N-1):
        nTe = i * Te
        signal_sum += sampled_signal[i] * sinc(t-nTe)
    return signal_sum

shannon = np.vectorize(shannon_reconstruct)

x = np.linspace(0, 2 * np.pi, N-1)
x_analog = np.linspace(0, 2*np.pi, 1000)

pl.plot(x0, sampled_signal, 'r')
pl.plot(x_analog, shannon(x_analog), 'y')
pl.figure()
pl.plot(x_analog, mySignal(x_analog) - shannon(x_analog))
pl.show()
../../_images/TP2 - FFT - Exercices_8_0.png ../../_images/TP2 - FFT - Exercices_8_1.png

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

import imageio as imio
colored_image = imio.imread('Moire.jpg')
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)
(599, 493, 3)
<matplotlib.image.AxesImage at 0x7f59e8f79a90>
../../_images/TP2 - FFT - Exercices_10_2.png
import imageio as imio
colored_image = imio.imread('services_0.jpg')
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)
(240, 365)
<matplotlib.image.AxesImage at 0x7f59e8e5a7d0>
<Figure size 432x288 with 0 Axes>
../../_images/TP2 - FFT - Exercices_11_3.png

2. Simple exercice with the spectrum of a signal.

A. On a 1D signal.

# Choose a grid x and a signal y
n = 1000
m = int(n/2)
print(m)
def step(x):
    return (x > 0) * 2 -1
x = np.linspace(-1,1,n)
# y = np.zeros_like(x)
# y[m:] = 1
# y[0:m] = -1
y = step(x)
plt.plot(x,y)
plt.title(" A step function ")
500
Text(0.5, 1.0, ' A step function ')
../../_images/TP2 - FFT - Exercices_14_2.png
# Use of FFT to compute the discrete Fourier transform
import scipy.fftpack as sf
spectre = sf.fft(y)
inverse_du_spectre = sf.ifft(spectre)

Ex: Compute the \(L^2\) of respectively the signal and its spectrum.

def normL2(y):
    yn = y/y.size # normalisation
    return np.sqrt(np.sum(np.abs(yn*yn)))

print('Norme L^2 du signal : ', normL2(y))
print('Norme L^2 du spectre :', normL2(spectre))
Norme L^2 du signal :  0.0316227766016838
Norme L^2 du spectre : 1.0000000000000002

Ex: Check numerically that the FFT and IFFT are numerical inverse of each others.

print('Norme de (inverse_du_spectre-spectre) : ', normL2(spectre-y), '<< 1, donc la FFT et IFFT sont bien inverse numérique')
Norme de (inverse_du_spectre-spectre) :  1.0004998750624612 << 1, donc la FFT et IFFT sont bien inverse numérique

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.

# illustration.
plt.figure(figsize = (20,8))
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. ")
Text(0.5, 1.0, ' The 0 frequency is in the middle. ')
../../_images/TP2 - FFT - Exercices_21_1.png
reduced = np.copy(spectre)
h = 40
reduced[h:n-h] = 0
reconstruct=sf.ifft(reduced)
plt.plot(np.real(reconstruct))
[<matplotlib.lines.Line2D at 0x7f59e8c5a610>]
../../_images/TP2 - FFT - Exercices_22_1.png

Ex: Read the script above and explain what it does. The illustrated phenomenon is called the Gibbs phenomenon and is due to the slow convergence of the Fourier serie at the discontinuites of the signal.

On obtient un signal filtré à l’aide d’un passe-bas (les hautes fréquences sont coupées).

plt.title("Decreasing of the spectrum log-modules")
fq = spectre[0:m] # + 00000000.1 # to avoid division by zero
plt.plot(np.log(np.abs(fq)))
/opt/hostedtoolcache/Python/3.7.9/x64/lib/python3.7/site-packages/ipykernel_launcher.py:3: RuntimeWarning: divide by zero encountered in log
  This is separate from the ipykernel package so we can avoid doing imports until
[<matplotlib.lines.Line2D at 0x7f59e8eaf590>]
../../_images/TP2 - FFT - Exercices_25_2.png

Ex: What can be observed on the shape of the spectrum ?

Plot odd and even coefficients. Can you explain why?

Why did we consider half of the signal?

odd = fq[1::2]
even = fq[0::2]

plt.figure(figsize = (20,8))
plt.subplot(1,2,1)
plt.plot(np.abs(odd))
plt.title(" Odd coeff ")
plt.subplot(1,2,2)
plt.plot(np.abs(even))
plt.title("Even coeff")
Text(0.5, 1.0, 'Even coeff')
../../_images/TP2 - FFT - Exercices_28_1.png

1 valeur sur 2 est nulle : c’est pour cela que l’on considère uniquement la moitié du signal

B. On an image.

import scipy.signal as sig
from scipy import misc
# 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)
<matplotlib.image.AxesImage at 0x7f59d355b250>
../../_images/TP2 - FFT - Exercices_32_1.png

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

spectre_im = sf.fft2(image_sub)
plt.figure(figsize = (20,8))
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))))
<matplotlib.image.AxesImage at 0x7f59d34b2190>
../../_images/TP2 - FFT - Exercices_34_1.png

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.

Les artefacts proviennent de la discontinuité lors du recollement de l’image (elle est sous échantillonnée).

def signal_porte(t):
    return ((t > -1/2) * (t < 1/2))

def gaussienne(x, sigma):
    return np.exp(-np.power(x, 2)/sigma**2)

tG = np.linspace(-1,1,100) #gauss
tD = np.linspace(-1,1,1000) # door

signal_gauss = gaussienne(tG, 0.5)
signal_gauss = signal_gauss / np.sum(signal_gauss) # normalization
porte = signal_porte(tD)
porte = porte / np.sum(porte) # normalization

def gaussianKernel(X,Y):
    K = np.zeros((X,Y))
    for x in range(X):
        for y in range(Y):
            x_ = x / (X-1) - 0.5
            y_ = y / (Y-1) - 0.5
            K[x,y] = gaussienne(x_,0.8) * gaussienne(y_, 0.8)
    return K / np.sum(K)

kernel = gaussianKernel(5,5)
plt.imshow(kernel)

denoised = sig.convolve2d(image_sub,kernel)
plt.figure(figsize = (20,8))
plt.subplot(1,2,1)
plt.gray()
plt.title(" Image denoise ")
plt.imshow(sig.convolve2d(image_sub,kernel))
plt.subplot(1,2,2)
plt.gray()
plt.title(" Image spectrum ")
plt.imshow(sf.fftshift(np.log(np.abs(sf.fft2(denoised)))))
<matplotlib.image.AxesImage at 0x7f59d33238d0>
../../_images/TP2 - FFT - Exercices_36_1.png ../../_images/TP2 - FFT - Exercices_36_2.png
# VERSION 2

def gaussDoor(t) :
    return np.convolve(gaussienne(np.linspace(-1, 1, 500), 0.8),signal_porte(t), 'same')

def makeKernel(X,Y, f):
    K = np.zeros((X,Y))
    for x in range(X):
        for y in range(Y):
            x_ = x / (X-1) - 0.5
            y_ = y / (Y-1) - 0.5
            K[x,y] = f[int(x_* len(f))] * f[int(y_ * len(f))]
    return K / np.sum(K)

# plt.plot(gaussDoor(np.linspace(-1, 1, 1000)))

kernel = makeKernel(5,5, gaussDoor(np.linspace(-1, 1, 1000)))
plt.imshow(kernel)

denoised = sig.convolve2d(image_sub,kernel)
plt.figure(figsize = (20,8))
plt.subplot(1,2,1)
plt.gray()
plt.title(" Image denoise ")
plt.imshow(sig.convolve2d(image_sub,kernel))
plt.subplot(1,2,2)
plt.gray()
plt.title(" Image spectrum ")
plt.imshow(sf.fftshift(np.log(np.abs(sf.fft2(denoised)))))
<matplotlib.image.AxesImage at 0x7f59d31e1050>
../../_images/TP2 - FFT - Exercices_37_1.png ../../_images/TP2 - FFT - Exercices_37_2.png

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.

import imageio as imio
bois = imio.imread('Moire.jpg')
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(sf.ifft2(new_spectrum)))
plt.subplot(1,3,3)
plt.title(" Larger random phase ")
plt.imshow(np.real(sf.ifft2(new_spectrum_2)))
<matplotlib.image.AxesImage at 0x7f59d3092110>
../../_images/TP2 - FFT - Exercices_39_1.png

Ex: Redo the experiment with the animal image. What do you observe?

racoon = misc.face()
texture_im = np.sum(racoon*[ 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.figure(figsize = (20,8))
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(sf.ifft2(new_spectrum)))
plt.subplot(1,3,3)
plt.title(" Larger random phase ")
plt.imshow(np.real(sf.ifft2(new_spectrum_2)))
<matplotlib.image.AxesImage at 0x7f59d2f3b890>
<Figure size 720x720 with 0 Axes>
../../_images/TP2 - FFT - Exercices_41_2.png

Le raton laveur est flouté. On a un bruit qui a été appliqué.

3. Convolution.

Ex: Construct a gaussian filter of size \((256,256)\) and make a convolution of the animal image with this filter. Then, show the image and zoom in a particular region to better see the effect of the convolution when the filter has a small standard deviation.

### Construct a gaussian filter.

def gaussienne(x, sigma):
    return np.exp(-np.power(x, 2)/sigma**2)

t = np.linspace(-1,1,100)
signal_gauss = (gaussienne(t, 0.5))
signal_gauss = signal_gauss / np.sum(signal_gauss)

def gaussianKernel(X,Y):
    K = np.zeros((X,Y))
    for x in range(X):
        for y in range(Y):
            x_ = x / (X-1) - 0.5
            y_ = y / (Y-1) - 0.5
            K[x,y] = gaussienne(x_,20) * gaussienne(y_, 20)
    return K / np.sum(K)

kernel = gaussianKernel(256,256)
fig = plt.figure(figsize=(15,15))
plt.imshow(kernel)
<matplotlib.image.AxesImage at 0x7f59d2e88790>
../../_images/TP2 - FFT - Exercices_45_1.png
fig = plt.figure(figsize=(15,15))
plt.gray()  # show the filtered result in grayscale
plt.subplot(1,2,1)
plt.title(" Original image ")
plt.imshow(image)

convoluted_image = sig.convolve2d(image,kernel)

plt.subplot(1,2,2)
plt.title(" Convoluted image ")
plt.imshow(convoluted_image)
plt.show()
plt.figure(figsize=(15,15))
plt.subplot(1,2,2)
plt.title(" Convoluted image ")
plt.imshow(convoluted_image[:200, :200])
plt.show()
../../_images/TP2 - FFT - Exercices_47_0.png

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.

from scipy.ndimage import gaussian_filter

sigma = 1
arch = imio.imread('services_0.jpg')

fig = plt.figure(figsize=(15,15))
plt.subplot(1,2,1)
plt.title("Original image")
plt.imshow(arch)

convoluted = gaussian_filter(arch, sigma=sigma)
sub_sampled = convoluted[::2,::2]

plt.subplot(1,2,2)
plt.title("Sub sampled image")
plt.imshow(sub_sampled)
plt.show()
../../_images/TP2 - FFT - Exercices_49_0.png

3. Linear and non-linear approximation.

### 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 ratio (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)
#def LinearApproximation(x,n):
def lowPassFrequency(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 = lowPassFrequency(image_sub,60)
approx2 = lowPassFrequency(image_sub,115)
plt.figure(figsize = (15,15))
plt.subplot(1,3,1)
plt.title("psnr: " + str(round(psnr(image_sub,approx1), 3)) + ", snr: " + str(round(snr(image_sub,approx1), 3)))
plt.imshow(approx1)
plt.subplot(1,3,2)
plt.title("psnr: " + str(round(psnr(image_sub,approx2), 3)) + ", snr: " + str(round(snr(image_sub,approx2), 3)))
plt.imshow(approx2)
plt.subplot(1,3,3)
plt.title(" Initial image ")
plt.imshow(image_sub)
<matplotlib.image.AxesImage at 0x7f1d15f09d00>
../../_images/TP2 - FFT - Exercices_52_1.png

Ex. Explain what is done in the previous cell.

Une image est filtrée dans le domaine fréquenciel (filtre passe bas). On perd donc de l’information dans l’image. Vient ensuite le calcul du ration signal/bruit pour quantifier la perte d’information.

Ex: Redo the experiment but in considering the non-linear approximation of the image by taking the first \(M\) larger coefficients in the spectrum.

Compare the linear and the non-linear approximation.

def NonLinearApproximation(x, m):
    spectre = sf.fftshift(sf.fft2(x))

    # plt.imshow(sf.fftshift(np.log(np.abs(spectre))))
    #fin Mth larger coefficient
    flatten = spectre.flatten()
    flatten.sort()
    flatten = flatten[::-1]
    mth = flatten[m]
    filteredSpectre = np.copy(spectre)
    
#     for i in range(filteredSpectre.shape[0]):
#         for j in range(filteredSpectre.shape[1]):
#             if(filteredSpectre[i, j] < mth):
#                 filteredSpectre[i, j] = 0
    # more concise
    filteredSpectre = np.multiply((filteredSpectre > mth), filteredSpectre)

    return np.real(sf.ifft2(sf.fftshift(filteredSpectre)))
approx1 = NonLinearApproximation(image_sub, int(image_sub.size*1)-1)
approx2 = NonLinearApproximation(image_sub, 4)
plt.figure(figsize = (15,15))
plt.subplot(1,3,1)
plt.title("psnr: " + str(round(psnr(image_sub,approx1), 3)) + ", snr: " + str(round(snr(image_sub,approx1), 3)))
plt.imshow(approx1)
plt.subplot(1,3,2)
plt.title("psnr: " + str(round(psnr(image_sub,approx2), 3)) + ", snr: " + str(round(snr(image_sub,approx2), 3)))
plt.imshow(approx2)
plt.subplot(1,3,3)
plt.title(" Initial image ")
plt.imshow(image_sub)
<matplotlib.image.AxesImage at 0x7f1d16140df0>
../../_images/TP2 - FFT - Exercices_56_1.png

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.

# utils to load the sounds.
import numpy as np
import wave as wv

def load_sound(file,n0):
    
    from scipy.io.wavfile import read
    s, data = read(file)
    x = np.array(data, dtype=float)
    x = x[:, 1] # keep only mono
    #print(x.shape)

#     x_raw = wv.open(file)
#     n = x_raw.getnframes()
#     x = np.frombuffer(x_raw.readframes(-1), 'Int16')
#     x_raw.close()

    if n0 !=0 and n0 < x.size:
        x = x[:n0]
    
    print(x)
    return x/np.max(x)

Hereafter, we load the \(3\) sounds and plot the second one.

n = 1024*16
s = 3 #number of signals.
x = np.zeros([n,3])
x[:,0] = load_sound("bird.wav", n)
x[:,1] = load_sound("female.wav", n)
# x[:,2] = load_sound("male.wav", n)

plt.figure(figsize = (10,10))
plt.plot(x[:, 1])
[-6.000e+00 -4.000e+00 -3.000e+00 ...  3.173e+03  2.311e+03  8.280e+02]
[-112. -229. -292. ... -713. -679. -676.]
[<matplotlib.lines.Line2D at 0x7f1d15e67dc0>]
../../_images/TP2 - FFT - Exercices_61_2.png

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

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.

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)))
<matplotlib.image.AxesImage at 0x7f1d1575e640>
../../_images/TP2 - FFT - Exercices_65_1.png

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

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))
7.663305028034717e-29

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.

nbre_selec = 600
from random import shuffle
print(np.shape(y)[0])
liste = list(range(np.shape(y)[0]))
shuffle(liste)
plt.plot(y[liste[0:nbre_selec],0],y[liste[0:nbre_selec],1],"o",markersize=1)
16384
[<matplotlib.lines.Line2D at 0x7f1d156769a0>]
../../_images/TP2 - FFT - Exercices_69_2.png

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

nbre_selec = 3000
H = np.asarray([np.imag(w1.flatten()),np.imag(w2.flatten())]).transpose()
liste2 = list(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)
[<matplotlib.lines.Line2D at 0x7f1d155d4520>]
../../_images/TP2 - FFT - Exercices_71_1.png

Les points ont l’air d’être réparti le long de deux “droites”.

Q. Comment the code below after running it.

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))
(16705,)
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()
../../_images/TP2 - FFT - Exercices_75_0.png
## Affiche les valeurs d'angles les plus importantes.
theta = []
for xx in (h>0.01)*t:
    if xx>0:
        theta.append(xx)
print(theta)
[0.1802282101269934, 0.1880625075964717, 0.19589680506595, 0.20373110253542828, 1.018498039361172, 1.0263323368306503, 1.0341666343001286, 1.0420009317696068, 1.0498352292390851, 1.0576695267085636, 1.065503824178042]

Q. Comment the code below after running it ?

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)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-40-7469da7aea16> in <module>
----> 1 projections = np.dot(Estimated,W)
      2 C = np.abs(projections)
      3 temp = np.max(C,0)
      4 I = np.argmax(C,0)
      5 threshold = .005

NameError: name 'Estimated' is not defined

La variable Estimated n’existe pas. Cela pose problème pour continuer.

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))
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,:])))
print(np.shape(MaListe[1]))
i = 1
plt.figure(figsize = (20,20))
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])))
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))
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]))