pulseanalyzer.cpp

/***************************************************************************
        Copyright (C) 2002-2015 Kentaro Kitagawa
                           kitagawa@phys.s.u-tokyo.ac.jp

        This program is free software; you can redistribute it and/or
        modify it under the terms of the GNU Library General Public
        License as published by the Free Software Foundation; either
        version 2 of the License, or (at your option) any later version.

        You should have received a copy of the GNU Library General
        Public License and a list of authors along with this program;
        see the files COPYING and AUTHORS.
***************************************************************************/
//---------------------------------------------------------------------------
#include "pulseanalyzer.h"
#include "analyzer.h"

REGISTER_TYPE(XDriverList, NMRBuiltInNetworkAnalyzer, "NMR Built-In Network Analyzer");

XNMRBuiltInNetworkAnalyzer::XNMRBuiltInNetworkAnalyzer(const char *name, bool runtime,
        Transaction &tr_meas, const shared_ptr<XMeasure> &meas) :
        XSecondaryDriverInterface<XNetworkAnalyzer>(name, runtime, ref(tr_meas), meas),
        m_pulser(create<XItemNode<XDriverList, XPulser> >("Pulser", false, ref(tr_meas), meas->drivers(), true)),
        m_dso(create<XItemNode<XDriverList, XDSO> >("DSO", false, ref(tr_meas), meas->drivers(), true)),
        m_sg(create<XItemNode<XDriverList, XSG> >("SG", false, ref(tr_meas), meas->drivers(), true)) {
    connect(m_dso);
    connect(m_pulser);
    connect(m_sg);

    iterate_commit([=](Transaction &tr){
        const char *cand[] = {"OFF", "11", "21", "51", "101", "201", "401", "801", "1601", "3201", ""};
        for(const char **it = cand; strlen( *it); it++) {
            tr[ *points()].add( *it);
        }
        tr[ *this].m_sweeping = false;
    });
    XNetworkAnalyzer::start();
}
void
XNMRBuiltInNetworkAnalyzer::clear() {
    restart(CAL_NONE, true);
}
void
XNMRBuiltInNetworkAnalyzer::onCalOpenTouched(const Snapshot &shot, XTouchableNode *) {
    restart(CAL_OPEN);
}
void
XNMRBuiltInNetworkAnalyzer::onCalShortTouched(const Snapshot &shot, XTouchableNode *) {
    restart(CAL_SHORT);
}
void
XNMRBuiltInNetworkAnalyzer::onCalTermTouched(const Snapshot &shot, XTouchableNode *) {
    restart(CAL_TERM);
}
void
XNMRBuiltInNetworkAnalyzer::onCalThruTouched(const Snapshot &shot, XTouchableNode *) {
    restart(CAL_THRU);
}
void
XNMRBuiltInNetworkAnalyzer::onStartFreqChanged(const Snapshot &shot, XValueNodeBase *) {
    clear();
}
void
XNMRBuiltInNetworkAnalyzer::onStopFreqChanged(const Snapshot &shot, XValueNodeBase *) {
    clear();
}
void
XNMRBuiltInNetworkAnalyzer::onAverageChanged(const Snapshot &shot, XValueNodeBase *) {
    clear();
}
void
XNMRBuiltInNetworkAnalyzer::onPointsChanged(const Snapshot &shot, XValueNodeBase *) {
    clear();
    int pts = atoi(Snapshot( *this)[ *points()].to_str().c_str());
    if( !pts){
        try {
            startContSweep();
        }
        catch (XInterface::XInterfaceError &e) {
            gErrPrint(e.msg());
        }
//      stop();
    }
}
void
XNMRBuiltInNetworkAnalyzer::getMarkerPos(unsigned int num, double &x, double &y) {
    Snapshot shot( *this);
    switch(num) {
    case 0:
        x = shot[ *this].m_marker_min.first;
        y = shot[ *this].m_marker_min.second;
        break;
    case 1:
        x = shot[ *this].m_marker_max.first;
        y = shot[ *this].m_marker_max.second;
        break;
    default:
        throw XDriver::XSkippedRecordError(__FILE__, __LINE__);
    }
}
void
XNMRBuiltInNetworkAnalyzer::oneSweep() {
    bool ret = restart(CAL_NONE);
    if( !ret)
        throw XDriver::XSkippedRecordError(__FILE__, __LINE__);
    while(Snapshot( *this)[ *this].m_sweeping) {
        msecsleep(30);
    }
}
bool
XNMRBuiltInNetworkAnalyzer::restart(int calmode, bool clear) {
    bool ret = false;
    iterate_commit([=, &ret](Transaction &tr){
        try {
            ret = false;
            restart(tr, calmode, clear);
            ret = true;
        }
        catch (XDriver::XSkippedRecordError &) {
        }
        catch (XInterface::XInterfaceError &e) {
            gErrPrint(e.msg());
        }
    });
    return ret;
}
void
XNMRBuiltInNetworkAnalyzer::restart(Transaction &tr, int calmode, bool clear) {
    Snapshot &shot_this(tr);

    int pts = atoi(shot_this[ *points()].to_str().c_str());
    tr[ *this].m_sweepPoints = pts;
    if( !pts)
        throw XDriver::XSkippedRecordError(__FILE__, __LINE__);

    tr[ *this].m_ftsum.clear();
    tr[ *this].m_ftsum_weight.clear();
    tr[ *this].m_calMode = calmode;
    if(clear || tr[ *this].m_raw_open.empty()) {
        tr[ *this].m_raw_open.clear();
        tr[ *this].m_raw_short.clear();
        tr[ *this].m_raw_term.clear();
        tr[ *this].m_raw_thru.clear();
        tr[ *this].m_raw_open.resize(pts, 1.0);
        tr[ *this].m_raw_short.resize(pts, -1.0);
        tr[ *this].m_raw_term.resize(pts, 0.0);
        tr[ *this].m_raw_thru.resize(pts, 1.0);
    }

    shared_ptr<XPulser> pulse = shot_this[ *m_pulser];
    if( !pulse)
        throw XDriver::XSkippedRecordError(__FILE__, __LINE__);
    shared_ptr<XSG> sg = shot_this[ *m_sg];
    if( !sg)
        throw XDriver::XSkippedRecordError(__FILE__, __LINE__);
    shared_ptr<XDSO> dso = shot_this[ *m_dso];
    if( !dso)
        throw XDriver::XSkippedRecordError(__FILE__, __LINE__);

    Snapshot shot_dso( *dso);
    double interval;
    int dso_len;
    if( !shot_dso[ *dso].time() || !shot_dso[ *dso].numChannels()) {
        interval = 1e-6; //temporary
        dso_len = 10000; //temporary
    }
    else {
        interval = shot_dso[ *dso].timeInterval();
        dso_len = shot_dso[ *dso].length();
    }

    double fmax = shot_this[ *stopFreq()];
    tr[ *this].m_sweepStop = fmax;
    double fmin = shot_this[ *startFreq()];
    tr[ *this].m_sweepStart = fmin;
    if((fmax <= fmin) || (fmin <= 0.1))
        throw XDriver::XSkippedRecordError(i18n("Invalid frequency settings"), __FILE__, __LINE__);

    pulse->iterate_commit([=, &tr](Transaction &trp){
        double plsbw = trp[ *pulse].paPulseBW() * 1e-3; //[MHz]
        double fstep = plsbw *1.0;
        fstep = std::max(fstep, (fmax - fmin) / (pts - 1));
        if(fstep < 0.001)
            throw XDriver::XSkippedRecordError(i18n("Invalid frequency settings"), __FILE__, __LINE__);
        tr[ *this].m_sweepStep = fstep;
        trp[ *pulse->pulseAnalyzerMode()] = true;
        double rept_ms = std::max(2.0 / ((fmax - fmin) / (pts - 1)) * 1e-3 * 2, 0.2);
        rept_ms = interval * 1e3 * lrint(rept_ms / (interval * 1e3)); //round to DSO interval.
        trp[ *pulse->paPulseRept()] = rept_ms;
        trp[ *pulse->output()] = true;
    });

    trans( *sg->freq()) = fmin;

    dso->iterate_commit([=](Transaction &trd){
        int avg = std::max(1L, lrint(0.03 / (interval * dso_len)));
        avg *= std::max(1u, (unsigned int)shot_this[ *average()]);
        trd[ *dso->average()] = (avg + 3) / 4 * 4; //round to phase cycling for NMR.
        trd[ *dso->firEnabled()] = false;
        trd[ *dso->restart()].touch(); //Restart averaging in DSO.
    });

    tr[ *this].m_sweeping = true;
}
void
XNMRBuiltInNetworkAnalyzer::startContSweep() {
    Snapshot shot_this( *this);
    shared_ptr<XPulser> pulse = shot_this[ *m_pulser];
    if(pulse) {
        pulse->iterate_commit([=](Transaction &tr){
            tr[ *pulse->pulseAnalyzerMode()] = false;
            tr[ *pulse->output()] = false;
        });
    }
    shared_ptr<XSG> sg = shot_this[ *m_sg];
    if(sg) {
        trans( *sg->freq()) = (double)shot_this[ *this].m_marker_min.first;
    }
}
void
XNMRBuiltInNetworkAnalyzer::acquireTrace(shared_ptr<RawData> &, unsigned int ch) {

}
void
XNMRBuiltInNetworkAnalyzer::convertRaw(RawDataReader &reader, Transaction &tr) throw (XRecordError&) {
}
void
XNMRBuiltInNetworkAnalyzer::writeTraceAndMarkers(Transaction &tr) {
    Snapshot &shot_this(tr);
    double fmin = shot_this[ *this].m_sweepStart;
    double fmax = shot_this[ *this].m_sweepStop;
    int pts = shot_this[ *this].m_sweepPoints;
//  if( !pts)
//      throw XDriver::XSkippedRecordError(__FILE__, __LINE__);

    tr[ *this].m_startFreq = fmin;
    double df = (fmax - fmin) / (pts - 1);
    tr[ *this].m_freqInterval = df;
    tr[ *this].trace_().resize(pts);

    auto ftsum = &tr[ *this].m_ftsum[0];
    auto ftsum_weight = &tr[ *this].m_ftsum_weight[0];
    for(unsigned int i = 0; i < pts; i++) {
        ftsum[i] /= ftsum_weight[i];
    }

    //Stores calibration curves.
    switch(tr[ *this].m_calMode) {
    case CAL_NONE:
        break;
    case CAL_OPEN:
        tr[ *this].m_raw_open.resize(pts);
        std::copy(tr[ *this].m_ftsum.begin(), tr[ *this].m_ftsum.end(), tr[ *this].m_raw_open.begin());
        break;
    case CAL_SHORT:
        tr[ *this].m_raw_short.resize(pts);
        std::copy(tr[ *this].m_ftsum.begin(), tr[ *this].m_ftsum.end(), tr[ *this].m_raw_short.begin());
        break;
    case CAL_TERM:
        tr[ *this].m_raw_term.resize(pts);
        std::copy(tr[ *this].m_ftsum.begin(), tr[ *this].m_ftsum.end(), tr[ *this].m_raw_term.begin());
        break;
    case CAL_THRU:
        tr[ *this].m_raw_thru.resize(pts);
        std::copy(tr[ *this].m_ftsum.begin(), tr[ *this].m_ftsum.end(), tr[ *this].m_raw_thru.begin());
        break;
    }

    auto rawopen = &tr[ *this].m_raw_open[0];
    auto rawshort = &tr[ *this].m_raw_short[0];
    auto rawterm = &tr[ *this].m_raw_term[0];
    auto trace = &tr[ *this].trace_()[0];
    for(unsigned int i = 0; i < pts; i++) {
        auto zport_in = - rawopen[i] / rawshort[i]; //Impedance of the port connected to LNA.
        auto s11_dut = 1.0 - 2.0 / ((1.0 + zport_in) / (1.0 - (ftsum[i] - rawterm[i]) / rawopen[i]) + 1.0 - zport_in);
        trace[i] = s11_dut;
    }

    //Tracking markers.
    auto &mkmin = tr[ *this].m_marker_min;
    mkmin.second = 1000;
    auto &mkmax = tr[ *this].m_marker_max;
    mkmax.second = -1000;
    auto z = tr[ *this].trace();
    for(unsigned int i = 0; i < pts; i++) {
        double r = 20.0 * log10(std::abs( *z++));
        if(r < mkmin.second)
            mkmin = std::pair<double, double>(i * df + fmin, r);
        if(r > mkmax.second)
            mkmax = std::pair<double, double>(i * df + fmin, r);
    }
}
void
XNMRBuiltInNetworkAnalyzer::open() throw (XKameError &) {

}
bool
XNMRBuiltInNetworkAnalyzer::checkDependency(const Snapshot &shot_this,
    const Snapshot &shot_emitter, const Snapshot &shot_others,
    XDriver *emitter) const {
    if( !shot_this[ *this].m_sweeping)
        return false;
    shared_ptr<XPulser> pulse = shot_this[ *m_pulser];
    if( !pulse) return false;
    shared_ptr<XDSO> dso = shot_this[ *m_dso];
    if( !dso) return false;
    shared_ptr<XSG> sg = shot_this[ *m_sg];
    if( !sg) return false;
    if(emitter != dso.get())
        return false;
   if(shot_emitter[ *dso].timeAwared() < shot_others[ *sg].time()) return false;
    return true;
}
void
XNMRBuiltInNetworkAnalyzer::analyze(Transaction &tr, const Snapshot &shot_emitter, const Snapshot &shot_others,
    XDriver *emitter) throw (XRecordError&) {
    const Snapshot &shot_this(tr);
    int pts = shot_this[ *this].m_sweepPoints;
    if( !pts) {
        tr[ *this].m_sweeping = false;
        throw XDriver::XSkippedRecordError(__FILE__, __LINE__);
    }

    const Snapshot &shot_dso(shot_emitter);
    const Snapshot &shot_sg(shot_others);
    const Snapshot &shot_pulse(shot_others);
    shared_ptr<XPulser> pulse = shot_this[ *m_pulser];
    shared_ptr<XDSO> dso = shot_this[ *m_dso];
    shared_ptr<XSG> sg = shot_this[ *m_sg];

    assert(shot_dso[ *dso].time() );

    if(shot_dso[ *dso].numChannels() < 1) {
        throw XSkippedRecordError(i18n("No record in DSO"), __FILE__, __LINE__);
    }
    if(shot_dso[ *dso].numChannels() < 2) {
        throw XSkippedRecordError(i18n("Two channels needed in DSO"), __FILE__, __LINE__);
    }
    if( !shot_dso[ *dso->singleSequence()]) {
        g_statusPrinter->printWarning(i18n("Use sequential average in DSO."));
    }
    int dso_len = shot_dso[ *dso].length();

    double interval = shot_dso[ *dso].timeInterval();
    if (interval <= 0) {
        throw XSkippedRecordError(i18n("Invalid time interval in waveforms."), __FILE__, __LINE__);
    }
    int pos = lrint(shot_dso[ *dso].trigPos());
    if(pos >= dso_len) {
        throw XSkippedRecordError(i18n("Position beyond waveforms."), __FILE__, __LINE__);
    }
    if(pos < 0) {
        throw XSkippedRecordError(i18n("Position beyond waveforms."), __FILE__, __LINE__);
    }

    if(pulse) {
        if( !shot_pulse[ *pulse].isPulseAnalyzerMode())
            throw XSkippedRecordError(i18n("Pulser configured not in Built-In Network Analyzer Mode."), __FILE__, __LINE__);
    }
    double rept = shot_pulse[ *pulse].rtime() * 1e-3; //[s]
    double plsorg = shot_pulse[ *pulse].paPulseOrigin() * 1e-6; //[s]

    double fmin = shot_this[ *this].m_sweepStart;
    double fmax = shot_this[ *this].m_sweepStop;

    unsigned int fftlen = (unsigned int)floor(std::max(plsorg / interval, 1e-6 / ((fmax - fmin) / (pts - 1)) / interval * 2));
    fftlen = FFT::fitLength(std::min(fftlen, (unsigned int)floor(rept / interval)));
    if( !m_fft || m_fft->length() != fftlen) {
        m_fft.reset(new FFT(-1, fftlen));
        m_fftin.resize(fftlen);
        m_fftout.resize(fftlen);
    }
    std::fill(m_fftin.begin(), m_fftin.end(), 0.0);
    unsigned int avg_in_wave = floor((dso_len - pos) / (rept / interval));
    if(avg_in_wave < 2)
        throw XSkippedRecordError(i18n("Too short waveforms."), __FILE__, __LINE__);
    avg_in_wave = avg_in_wave / 2 * 2;

    for(unsigned int av = 0; av < avg_in_wave; ++av) {
        int lpos = lrint(pos + (rept * av) / interval);
        int org = lrint(pos + (plsorg + rept * av) / interval);
        bool inverted = (av % 2 == 1);
        const double *wavecos = &shot_dso[ *dso].wave(0)[lpos];
        const double *wavesin = &shot_dso[ *dso].wave(1)[lpos];
        for(unsigned int i = fftlen - (org - lpos); i < fftlen; ++i) {
            m_fftin[i] += std::complex<double>( *wavecos++,  *wavesin++) * (inverted ? -1.0 : 1.0);
        }
        for(unsigned int i = 0; i < fftlen - (org - lpos); ++i) {
            m_fftin[i] += std::complex<double>( *wavecos++,  *wavesin++) * (inverted ? -1.0 : 1.0);
        }
    }
    double fft_df = 1.0 / (interval * fftlen) * 1e-6; //[MHz]
    m_fft->exec(m_fftin, m_fftout);

    double plsbw = shot_pulse[ *pulse].paPulseBW() * 1e-3; //[MHz]

    if(tr[ *this].m_ftsum.size() != pts) {
        tr[ *this].m_ftsum.resize(pts);
        tr[ *this].m_ftsum_weight.resize(pts);
        std::fill(tr[ *this].m_ftsum.begin(), tr[ *this].m_ftsum.end(), 0.0);
        std::fill(tr[ *this].m_ftsum_weight.begin(), tr[ *this].m_ftsum_weight.end(), 0);
    }
    double freq = shot_sg[ *sg].freq();
    tr[ *this].m_lastCenterFreq = freq;
    if(freq < fmin - plsbw/2) {
        restart(tr, shot_this[ *this].m_calMode);
        throw XDriver::XSkippedRecordError(__FILE__, __LINE__);
    }
    double normalize = 1.0 / avg_in_wave / fftlen;
    auto ftsum = &tr[ *this].m_ftsum[0];
    auto ftsum_weight = &tr[ *this].m_ftsum_weight[0];
    for(int i = 0; i < fftlen; ++i) {
        double f = fft_df * ((i >= fftlen / 2) ? i - (int)fftlen : i);
        if(abs(f) > plsbw / 2)
            continue;
        f += freq;
        int j = lrint((f - fmin) / (fmax - fmin) * (pts - 1));
        if((j < 0) || (j >= pts))
            continue;
        ftsum[j] += m_fftout[i] * normalize;
        ++ftsum_weight[j];
    }
    if( !shot_this[ *this].m_ftsum.size()) {
        restart(tr, shot_this[ *this].m_calMode);
        throw XDriver::XSkippedRecordError(__FILE__, __LINE__);
    }
    double fstep = shot_this[ *this].m_sweepStep;
    if(freq + fstep / 2 > fmax) {
        writeTraceAndMarkers(tr);
        tr[ *this].m_sweeping = false;
    }
    throw XDriver::XSkippedRecordError(__FILE__, __LINE__);
}
void
XNMRBuiltInNetworkAnalyzer::visualize(const Snapshot &shot) {
    if(shot[ *this].m_sweeping) {
        double freq = shot[ *this].m_lastCenterFreq;
        double fstep = shot[ *this].m_sweepStep;

        freq += fstep;
        shared_ptr<XSG> sg = shot[ *m_sg];
        assert(sg);
        trans( *sg->freq()) = freq; //setting new freq.
    }
    else
        XNetworkAnalyzer::visualize(shot);
}