CODEWITHSUNDEEP
pythonnumpyscipy
samane filabi
samane filabi

Reputation: 15

Phase Shift between two noisy signals of Stochastic Resonance

for task 3 I wanna try to calculate phase shift between my noisy signal and Periodic forcing, I used Cross correlation method but phase shift becomes very small and does not make sense, I filtered my noisy signal or different methods but did not work , do you have any Idea how I should calculate exactly phase shift ?Image of Task

At first I plotted Noisy signal and than Calculated phase shift between noisy signal(x) and periodic forcing (x2)


import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as signal
simulation_time = 100.0
sigma = 0.3
A =10
T=100
omega = 0.1
dt = 0.001
sqdt = np.sqrt(dt) #Precompute
time = np.arange(0, simulation_time, dt)
x = np.empty(len(time))
x2 = np.empty(len(time))
#Initial conditions 
x[0] = 1.0

i = 1
for t in time[1:]:
    fx = x[i-1] - x[i-1]*x[i-1]*x[i-1] + A*np.cos(omega*t)
    x2[i]=A*np.cos(omega*t)
    gx = sigma
    x[i] = x[i-1] + dt*fx + sqdt*gx*np.random.standard_normal()
    i += 1
plt.figure()
plt.plot(time, x)
plt.plot(time, x2)
plt.xlabel('time')
plt.ylabel('signal')
plt.show()

#cross corelation method 
output = np.correlate(x,x2*0.15,mode='same')
lags   = time-50
plt.plot(lags,output)
maximum = np.max(output)
phase_shift = lags[output==maximum][0]
phase_shift


0.05799999999999983

Upvotes: 0

Views: 154

Answers (1)

Reinderien
Reinderien

Reputation: 15273

When in doubt, verify against an FFT. This is not as efficient, but is more intuitive to interpret:

import numpy as np
import matplotlib.pyplot as plt


simulation_time = 100.0
sigma = 0.3
A = 10
T = 100
omega = 0.1
dt = 0.001
sqdt = np.sqrt(dt)  # Precompute
time = np.arange(0, simulation_time, dt)
x = np.empty(len(time))
x2 = np.empty(len(time))

# Initial conditions
x[0] = 1.0

i = 1
for t in time[1:]:
    fx = x[i-1] - x[i-1]*x[i-1]*x[i-1] + A*np.cos(omega*t)
    x2[i] = A*np.cos(omega*t)
    gx = sigma
    x[i] = x[i-1] + dt*fx + sqdt*gx*np.random.standard_normal()
    i += 1


spectrum = np.fft.rfft(x)
freqs = np.fft.rfftfreq(len(x), d=dt)
i_peak = np.argmax(spectrum)
phase = np.angle(spectrum[i_peak])
shift = phase+np.pi/2

print(f'Peak frequency: {freqs[i_peak]} Hz - should be close to {omega/2/np.pi:.6}')
print(f'Peak phase: {phase:.6} rad')
print(f'Phase shift from cosine: {shift:.3} rad, or {shift/2/np.pi:.1%}')

plt.figure()
plt.plot(time, x)
plt.plot(time, x2)
plt.xlabel('time')
plt.ylabel('signal')
plt.show()

Peak frequency: 0.02 Hz - should be close to 0.0159155
Peak phase: -1.4046 rad
Phase shift from cosine: 0.166 rad, or 2.6%

Upvotes: 2

Related Questions