Spectral Analysis

Auto-correlation Function

Auto-correlation Function

Lets take the following signal

s(t) = \frac{1}{1 - 0.9 sin (t)}

Question

  • Compute the discrete Fourier transform of the signal

  • Compute and plot in log-log the power spectral

  • Add some random noise to the signal using a random number generator: . Where is an adjustable parameter.

  • Plot the noisy signal and the corresponding PSD.

  • Compare the PSD of the original signal and the noisy one (at small scale)

Hint

To compute the PSD use the fft function implemented in numpy module (python)

Solution
1
"""
2
@author: yacine.mezemate
3
"""
4
import scipy as sp
5
import numpy as np
6
import matplotlib.pylab as plt
7
import random
8
import math
9
10
max = 4000
11
signal = sp.zeros((max),float)
12
ps = sp.zeros((max),float)
13
step = 2*math.pi/1000
14
15
def func(array,max):
16
    f = sp.zeros((max + 1),float) 
17
    step = 2*sp.pi/1000; x = 0. 
18
    for i in range(0,max):
19
        f[i] = 1/(1. - 0.9*sp.sin(x))             # Function
20
        signal[i] = (1./(1-0.9*sp.sin(x)))+0.5*(2.*random.random()-1.) # noise
21
        x += step
22
        
23
def filtr():                           # Low-pass windowed - sinc filter
24
    y = sp.zeros((max),float); h = sp.zeros((max),float)
25
    m = 100                              # Set filter length (101 points)
26
    fc = .07   
27
    for i in range(0,100):             # Calculate low-pass filter kernel
28
        if ((i-(m/2)) == 0):  h[i] = 2*sp.pi*fc 
29
        if ((i-(m/2))!= 0):   h[i] = sp.sin(2*sp.pi*fc*(i-m/2))/(i-m/2)
30
        h[i] = h[i]*(0.54 - 0.46*sp.cos(2*sp.pi*i/m))          # Hamming window
31
    sum = 0.                           # Normalize low-pass filter kernel
32
    for i in range(0,100):  sum = sum + h[i]
33
    for i in range(0,100):  h[i] = h[i] / sum 
34
    for j in range(100,max-1):               # Convolve input with filter
35
        y[j] = 0.                    
36
        for i in range(0,100): y[j] = y[j] + signal[j-i] * h[i]
37
    
38
    return y
39
    
40
#Calulate the PSD
41
def spectra(data):
42
    
43
    L = sp.power(2,12)
44
    data = data[0:L]
45
    FMax=np.floor(L/2)
46
    f = sp.arange(1,FMax+1,1)
47
    E = sp.zeros((FMax,))
48
    
49
    Buffer = np.absolute(np.fft.fft(data))
50
    E = Buffer[0:FMax]*Buffer[0:FMax]
51
    plt.figure(0)
52
    plt.loglog(f,E)
53
    plt.xlabel("f", fontsize =18)
54
    plt.ylabel("PSD", fontsize =18)
55
        
56
func(signal,max)
57
spectra(signal)
58
59
fi = filtr()
60
spectra(fi)
61
62
63
#plot noise signal
64
plt.figure(1)
65
plt.plot(signal)
66
plt.xlabel("t", fontsize =18)
67
plt.ylabel("y(t)", fontsize =18)
68
69
#plot filtred signal
70
plt.figure(2)
71
plt.plot(fi)
72
plt.xlabel("t", fontsize =18)
73
plt.ylabel("s(t)", fontsize =18)
(a) Noisy signal, (b): Denoised signal, (c): The corresponding PSD
PreviousPreviousNextNext
HomepageHomepagePrintPrintCreated with Scenari (new window)