%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import constants
from scipy import signal as sigproc
from scipy.fft import next_fast_len

shot_no = 36649

Data loading

The data is read from ISF files, using the isfreader library. It is installed bellow automatically if not already installed.

!pip3 install isfreader

Requirement already satisfied: isfreader in /opt/anaconda3/lib/python3.8/site-packages (1.0.0)
Requirement already satisfied: numpy>=1.9.1 in /opt/anaconda3/lib/python3.8/site-packages (from isfreader) (1.19.2)

import isfreader

ds = np.DataSource('/tmp')   # temporary sotarge for downloaded files
data_URL = 'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/Interferometry/DAS_raw_data_dir/ch{ch_id}.isf'

def load_channel(shot_no, chanel_id, channel_name=None, in_ms=True):
    fname = ds.open(data_URL.format(shot_no=shot_no, ch_id=chanel_id)).name
    data = isfreader.read_file(fname)
    signal = pd.Series(data[:,1], name=channel_name,
                       index=pd.Index(data[:,0], name='time'))
    if in_ms:
        signal.index *= 1e3  # s -> ms
    return signal

mixer = load_channel(shot_no, 1, 'mixer')

ref_saw = load_channel(shot_no, 3, 'ref_saw')

analog_dphase = load_channel(shot_no, 4, 'analog_dphase')

Sampling frequency in kHz

f_s = mixer.size / (mixer.index[-1] - mixer.index[0])  # ms -> kHz
f_s

250000.05000001

Spectral analysis of signals

The mixer signals is a base sine wave (the envelope of the mixing) at a frequency close to 500 kHz. The reference saw-tooth frequency sweeping wave has the same base frequency, but with a trail of harmonics forming the sharp saw-tooth shape.

def show_spectrum(signal, target_dfreq=10):
    nperseg = int(f_s / target_dfreq) 
    f, psd = sigproc.welch(signal.values, fs=f_s,
                           nperseg=nperseg, nfft=next_fast_len(nperseg))
    
    f_max = f[psd.argmax()]
    ax = plt.gca()
    ax.plot(f, psd, label=signal.name)
    ax.set(xlabel='frequency [kHz]', ylabel='power spectral density [V^2]')
    return f_max

show_spectrum(ref_saw)
f_base = show_spectrum(mixer)
plt.axvline(f_base, label=f'base frequency f={f_base:.0f} kHz', color='C3')
plt.loglog()
plt.legend()
f_base

540.0001080000216

Base band signal analysis

The difference between the phase of the reference and mixer signal is proportional to the line-integrated electron density. Therefore, it is necessary to obtain the "phase" of both signals.

The instantanous phase and amplitude of the base signal can be inferred only for the base bad, i.e. by removing higher and lower frequencies (i.e. extracting the base sine wave from the reference saw-tooth signal)

base_band_hwidth = 50 # kHz

base_band_filter = sigproc.iirfilter(8, 
                    [f_base - base_band_hwidth, f_base+base_band_hwidth],
                                     fs=f_s,
                                    btype='bandpass', output='sos')

def freq_filter(signal, sos_filter):
    return pd.Series(name=signal.name+'_filtered', index=signal.index,
        data=sigproc.sosfiltfilt(sos_filter, signal.values))

mixer_filtered = freq_filter(mixer, base_band_filter)
ref_saw_filtered = freq_filter(ref_saw, base_band_filter)

The instantanous phase and amplitude of this base band signal can be obtained from the analytic signal using the Hilbert transform. For convenience of processing and visualization, the phase is unwrapped ($\pi$ discontinuities are shifted) and it's time derivative, the instantaneous frequency, is calucated as well.

def inst_phase_amp(signal):
    analytic_signal = sigproc.hilbert(signal.values)
    inst_phase = pd.Series(np.unwrap(np.angle(analytic_signal)), 
                      index=signal.index, 
                     name=signal.name+'_inst_phase')
    # remove base band phase evolution trend
    inst_freq = pd.Series(np.gradient(inst_phase, inst_phase.index), 
                          index=inst_phase.index, name=signal.name+'_inst_freq')
    inst_amp = pd.Series(np.abs(analytic_signal), index=signal.index, 
                     name=signal.name+'_inst_amp')
    return inst_phase, inst_freq, inst_amp

mixer_phase, mixer_freq, mixer_amp = inst_phase_amp(mixer_filtered)
ref_saw_phase, ref_saw_freq, ref_saw_amp = inst_phase_amp(ref_saw_filtered)

For convenience of visualization it is convenient to subtract the phase evolution of the base nad signal from the phases (which is just a linear trend). This base trend vanishes in the final phase difference calculation anyways as long as the same trend is removed. Therefore, the referense saw-tooth signals is used for the trend estimation.

omega_base_ref = ref_saw_freq.median()
mixer_dfreq = mixer_freq - omega_base_ref
mixer_dphase = mixer_dfreq.cumsum() / f_s
ref_saw_dphase = (ref_saw_freq - omega_base_ref).cumsum() / f_s
omega_base_ref / (2*np.pi)

539.0220549190738

mixer_freq.loc[-1:0].median() / (2*np.pi)

539.0027640932642

The signal boundaries may be distorted by the filter, so they are not used.

boundary_Dt = 0.5  # in ms
in_boundary = slice(mixer_dfreq.index[0]+boundary_Dt, 
                    mixer_dfreq.index[-1] - boundary_Dt)
in_boundary

slice(-1.5019999999999998, 17.497996, None)

ref_saw_dphase.loc[in_boundary].plot()
mixer_dphase.loc[in_boundary].plot()
plt.legend();

<matplotlib.legend.Legend at 0x7faf6b37ce20>

The phase difference between the interferometer signal with respect to the reference signal is proportional to the line-averaged density. An arbitrary initial phase difference is subtracted.

dphase = pd.Series(ref_saw_dphase - mixer_dphase, 
                   index=mixer.index, name='dphase')
dphase -= dphase.loc[-1:0].mean()

analog_dphase -= analog_dphase.loc[-1:0].mean()
analog_dphase *= 3*np.pi # analog conversion?

fig, axs = plt.subplots(2, 1, sharex=True, figsize=(8, 8))
plt.sca(axs[0])
mixer_amp.loc[in_boundary].plot()
axs[0].set(ylabel='mixer amplitude', title=f'GOLEM #{shot_no}')
plt.sca(axs[1])
analog_dphase.plot()
dphase.loc[in_boundary].plot()
plt.legend()
plt.savefig('icon-fig.png')