# Simulated example¶

This is a clean simulated deconvolution example. That is, using the measured frequency response an output signal is generated by means of convolution and thereafter the input signal estimated by deconvolution. The advantage of this approach is to verify the general applicability of the method.

In [1]:
%pylab inline
from GUM2DFT import GUM_DFT,GUM_iDFT,GUMdeconv, AmpPhase2DFT
from statsmodels.tsa.arima_process import ArmaProcess

Populating the interactive namespace from numpy and matplotlib


### Measurement system¶

In [9]:
calib = np.loadtxt("data/calib_PA055.dat")
f = calib[:,0]
Nf = 2*len(f)-1
FR = calib[:,1]*np.exp(1j*calib[:,3])
uAmp = calib[:,2]
uPhas= calib[:,4]


### Input signal and simulated output signal¶

In [3]:
input_file = "data/reference_signal.dat"
noise_std = 1e-4

Ts = 2e-9

P = np.fft.rfft(np.r_[pressure,np.zeros(2*len(f)-len(pressure)-1)])
x = np.fft.irfft(FR*P)[:len(time)]

In [4]:
# adding correlated noise (ARMA process)
N = len(time)
ar = np.array([0.75,-0.6])
ma = np.array([0.45,0.1])
ARMA = ArmaProcess(ar,ma)
noise = ARMA.generate_sample(N)*0.01

x += np.random.randn(len(time))*noise_std

acov = ARMA.acovf(nobs=N)*noise_std**2

Ux = np.zeros((N,N))
for k in range(N):
Ux[k,k:] = acov[:N-k]
Ux = Ux + Ux.T - np.diag(np.diag(Ux))

In [8]:
figure(figsize=(14,8));clf()
plot(time*1e6,pressure,label="true pressure signal")
plot(time*1e6,x,label="simulated measurement")
legend(loc="best");
xlabel(r"time / $\mu$s",fontsize=20)
ylabel("pressure / MPa",fontsize=20)

Out[8]:
<matplotlib.text.Text at 0x19b1d390>

### Regularized inverse system for deconvolution¶

In [10]:
# low-pass filter
fcut = 120e6
HL = 1/(1+1j*f/fcut)**2
uAmp /= np.abs(HL)
# inverse system for input estimation
H = FR/HL

In [12]:
figure(figsize=(14,12));clf()
subplot(211)
errorbar(f*1e-6,calib[:,1],2*uAmp,fmt=".",alpha=0.2)
xlim(0.5,150);ylim(0,0.25)
xlabel("frequency / MHz",fontsize=22); tick_params(which="both",labelsize=18)
ylabel("amplitude / a.u.",fontsize=22)
subplot(212)
errorbar(f*1e-6,calib[:,2],2*uPhas,fmt=".",alpha=0.2)
xlim(0.5,150)
xlabel("frequency / MHz",fontsize=22); tick_params(which="both",labelsize=18)
ylabel("phase / deg",fontsize=22);


### Propagation of uncertainties¶

In [13]:
# propagation to real/imag domain
UH = AmpPhase2DFT(np.abs(FR),np.angle(FR)*np.pi/180,np.r_[uAmp,uPhas*np.pi/180])[1]

# propagation from time to frequency domain
Y,UY = GUM_DFT(x,Ux,N=Nf,return_complex=True)

# propagation through deconvolution equation
XH,UXH = GUMdeconv(H,Y,UH,UY)

# propagation back to time domain
xh,Uxh = GUM_iDFT(XH,UXH,Nx=len(x))

In [17]:
ux = np.sqrt(np.diag(Uxh))

figure(figsize=(14,8));clf()
plot(time*1e6,xh,label="estimate of measurand",linewidth=2)
plot(time*1e6,pressure,"--",label="true value of measurand",linewidth=2)
fill_between(time*1e6,xh+ux,xh-ux,alpha=0.2)
xlabel(r"time / $\mu$s",fontsize=22)
ylabel("signal amplitude / a.u.",fontsize=22)
tick_params(which="major",labelsize=18)
legend(loc="best",fontsize=20,fancybox=True)
xlim(0.5,1.5);

In [ ]: