displayroutines.cpp

USB 上位机程序 VNA使用,网络分析仪原理使用仪器

//
//    Copyright 2004-2007 Thomas C. McDermott, N5EG
//    This file is part of VNAR - the Vector Network Analyzer program.
//
//    VNAR is free software; you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation; either version 2 of the License, or
//    (at your option) any later version.
//
//    VNAR is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//    You should have received a copy of the GNU General Public License
//    along with VNAR, if not, write to the Free Software
//    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//

#pragma once
//
//
//	Routines to manipulate display points in rectangular
//  and polar format for the Vector Network Analyzer
//
//	7-4-03  TCMcDermott
//  7-10-03 Add Frequency_Grid Class.
//  7-11-03 Add CalDataSet.
//  7-12-03 Add routines to compute S-parameter error terms
//			  and convert measured S11 to actual S11.
//  7-13-03 Add routines to save CalDat to file, and load
//			  CalData from file.
//  7-27-03 Add routines to compute rectangular display
//			  coordinates.
//  8-17-03 Add Detector Class to support separate detector
//			  calibration.
//  9-16-03 Phase correction algorithm implemented for
//			  detector calibration.  ResolveTranPolar and
//			  ResolveReflPolar changed to use detector cal values.
//  3-28-04 Add TxLevel adjust to IDAC on the VNA board
//  4-15-04 Revise amplitude interpretation algorithm
// 11-08-04 Revise Detector amplitude calibration algorithm to account for
//			two sets of S21 readings - with and without 40 dB pad.
//			Also improve to one reading every dB of range for S21.
//			REFL detector does not have this done, but TRAN detector does,
//			so there are two different calibration algorithms. Also affects
//			MagTodB routine.  May have to zero out the LowMag routine for
//			REFL, but still take data every 1 dB like TRAN.
// 12-29-05 Add Coupler Ripple and Directivity calibration compensation
//			independent of fixture calibration. These are both controlled
//			by #define statements within this file. If they are changed, all
//			detector and fixture calibrations will need to be re-run.
// 03-20-06 Change MagTodBTran to use 3 measurements: Hi, Mid, and Lo
//
//



#include "stdafx.h"
#using <mscorlib.dll>
#using <System.dll>
using namespace System;
using namespace System::ComponentModel;
using namespace System::IO;
using namespace System::Windows::Forms;

#include <math.h>
#include <complex>
#include "DataDisplay.h"
#include "Constants.h"
#include "DisplayRoutines.h"

#define DIRECTCAL		// enable compensation of coupler directivity to magnitude
#define COUPLERCAL	// enable compensation of coupler ripple

	// resolve measured data sets to Magnitude and Phase

	// Draw reflection measurements on Polar Chart
	// Requires resolving the phase quandrant from I/Q measurements
	//   and magnitude from the reflection or tranmsission magnitude
	//	 (which is a logrithmic value).
	//
	// when ReflPI is low,  the phase is near   0 degrees
	// when ReflPI is high, the phase is near 180 degrees
	// when ReflPQ is low,  the phase is near  90 degrees
	// when ReflPQ is high, the phase is near -90 degrees
	//
	// The I and Q phase readings saturates near the bottom and top of the reading,
	//   but the other reading is then in the middle of it's useful range.
	//   This it's always possible to accurately resolve the phase and quadrant.
	//
	// Quadrant 1 is 0 to +90 degrees, quadrant 2 is +90 to +180,
	//   quadrant 3 is -180 to -90, quadrant 4 is 0 to -90 degrees.
	//

	/// Convert Mag and Phase to X,Y screen display coordinates in polar display mode
void ToDisplayPolar(double magnitude, double phase, int polarRad, int xoffset, int yoffset, int& X, int& Y)
{

	double fx, fy, phase_radians;
	int	x, y;

	phase_radians = phase / 180.0 * 3.1415926;

	fx = magnitude * cos(phase_radians);
	fy = magnitude * sin(phase_radians);

	// Scale to the size of the chart (polar radius)

	fx = fx * static_cast<double>(polarRad);
	fy = fy * static_cast<double>(polarRad);

	// Convert to integer

	x = (int)fx;
	y = (int)fy;

	// Center on the chart

	x = x + polarRad + xoffset;
	y = polarRad + yoffset - y;

	X = x;		// return X and Y display points
	Y = y;

}



	/// Convert Magnitude to rectangular (X, Y) display coordinates
int ToDisplayRectMag(double magnitude, int height, int dbScaleFactor, int refLevel)
{
	/// magnitude has a normalized value.
	/// Zero dB is the top of the screen - magnitude value = 1.000
	/// -100 dB is the bottom of the screen (@ 10 db/div) - magnitude value 1e-5
	/// Convert magnitude to vertical display accounting for the scale factor in dB/division

	double dbmag;

	if(magnitude < 0.0000001)
		dbmag = 140.0;
	else
		dbmag = refLevel - 20.0 * log10(magnitude);	// 0 to 100 for zero dB to -100 dB.

	dbmag /= ((double)dbScaleFactor * 10.0);			// 0 to +1 for 10db/div

	// scale to the vertical screen display area
	int vertical = (int)(dbmag * (double)height);	// 0 to scopeDisp.height for zero dB to height for max dB

	if (vertical > height)					// crop display - probably a better way to do this
		vertical = height;
	if (vertical < 0)
		vertical = 0;

	return(vertical);						// return Y as positive number

}
	/// Convert Phase to rectangular (X, Y) display coordinates
int ToDisplayRectPhs(double phase, int height)
{
	/// phase in degrees ranges from -180 to +180

	phase += 180.0;								///< offset the phase display so zero is centerline
	double range = phase * (double)height / 360.0; ///< top to bottom of display is 360 degrees

	return((int)range);
}
/// Convert value of groupdelay to Y-display coordinate
int ToDisplayRectGD(double groupdelay, int height, int scaleFactor)
{

	/// scale factor = 1  implies 100 picoseconds
	// convert to nanoseconds
	groupdelay *= 1000000000.0;

	/// set full screen display to 10 units of the selected scale
	// 100 picoseconds * 10 = 1 nanosecond full-scale at max resolution
	groupdelay /= scaleFactor;

	// scale it to display region
	groupdelay *= (double)height;

	/// clip display value to screen height
	if (groupdelay > height)
		groupdelay = height;

	return((int)groupdelay);
}
/// Construct a new frequency grid of size numPoints
FrequencyGrid::FrequencyGrid(int numPoints)
{
	if (numPoints > 1024)
		throw new System::IndexOutOfRangeException;

	FrequencyIndex = new int __gc[numPoints];
	startFreq = 200000;
	stopFreq = 120000000;
	points = numPoints;
	Build();
};
/// Set start frequency of grid
void FrequencyGrid::SetStartF(int start)
	{
		startFreq = start;
		Build();
	}
/// Set stop frequency of grid
void FrequencyGrid::SetStopF(int stop)
{
	stopFreq = stop;
	Build();
}
/// Convert gridpoint to it's frequency
int FrequencyGrid::Frequency(int gridpoint)
{
	/// if gridpoint is out of range, clip to actual frequency range

	if (gridpoint < 0)
		return (FrequencyIndex[0]);
	if (gridpoint >= points)
		return (FrequencyIndex[points-1]);

	return(FrequencyIndex[gridpoint]);
	
}
/// Derive DDS divisor value (as 64-bit integer) from Frequency
long long int FrequencyGrid::DDS(int Frequency)
{

/// Calculate the 48-bit synthesizer accumulator value from Frequency
/// the VNA crystal is ~12.000 MHz. It gets multiplied by 24 going to the DDS accumulator 
/// The clock rate is thus approximately 24 * ~12.000 MHz. = ~288.000 MHz.
/// The synthesizer accumulator value needed is:   N = Fdesired * 2^48 / Fclk

	// Modified -- N = Fdesired * 2^48 / Fclk.
	// Division done in two steps to prevent numeric overflow.
	// Calibration routine can determine a better value for Fclk than
	// our initial guess.
	
	// For calibration, emit F= 100.000000 Mhz.
	// Then new guess for Fclk = 2.88 * Fmeasured
	// Ferror is the Frequency calibration error at 100 MHz.
	// It has to be scaled up to 288 MHz.

	long long int N;
	N = Convert::ToInt64(Frequency) * 4294967296;  // Freq * 2^32
	N /= (288000000 + (Ferror * 288) / 100 );	   // Ferror * 2.88
	N *= 65536;								       // Freq * 2^16
	return(N);
}
/// Find nearest gridpoint to Frequency
/// \result is GridPoint Number or -1 if error
int FrequencyGrid::GridPoint(int Frequency)	
{

	if (((double)Frequency > stopFreq) || ((double)Frequency < startFreq))
		return(-1);				// Frequency is outside the grid range

	double result = ((double)Frequency - startFreq) / delta;
	int iresult = (int)(result + 0.5);

	return(iresult);
}
/* __property*/ int FrequencyGrid::get_Points() { return points; }		/// get number of points in grid
/* __property*/ int FrequencyGrid::get_StartF() { return startFreq; }	/// get start frequency of grid
/* __property*/ int FrequencyGrid::get_StopF() { return stopFreq; }		/// get stop frequency of grid
/* __property*/ int FrequencyGrid::get_Ferror() { return ferror; }		/// get Frequency error of the reference clock
/* __property*/ void FrequencyGrid::set_Ferror(int f) { ferror = f; }	/// set Frequency error of the reference clock
void FrequencyGrid::Build(void)
	{
		delta = ((double)stopFreq-(double)startFreq)/(points-1);

		// build the frequency grid
		for (int i = 0; i<points; i++)
			FrequencyIndex[i] = startFreq + (int)((double)i * delta);
	}

/// Detector Constructor, allocate default values
Detector::Detector(void)							
	{
		unitsPerdB = 57.0;				// default for uncalibrated data

		centerPhase = 1815;				// starting guess at phase center
		lowerPhase = 900;				// starting guess at lower and upper phase crossovers
		upperPhase = 2750;	
		centerQindex = 0, centerQ2index = 0;
		centerIindex = 0, centerI2index = 0, centerI3index = 0;
		quadratureError = 0.0;
		phaseSlope = 0;
		pherror = new double __gc[360];	// hold phaserror array
		MagTable = new int __gc[21, 4];	// hold raw refl Mag count data
		m = new double __gc[21];			// value of 'm' at each frequency
		b = new double __gc[21];			// value of 'b' at each frequency
		r = new double __gc[21];			// value of 'r' at each frequency
		flat = new double __gc[21];		// value of 'flat' at each frequency

		// New 10-17-2005
		DirMag = new int __gc[21];		// Directivity Magnitude ADC count at each frequnecy
		DirIphs = new int __gc[21];		// Directivity I phase ADC count at each frequency
		DirQphs = new int __gc[21];		// Directivity Q phase ADC count at each frequency

		calibrated = false;				// values are not yet calibrated
		DirCalibrated = false;			// directivity calibration not yet run

	}
	/// construct phase points from raw data
bool Detector::PhaseCal(int Iphase __gc[], int Qphase __gc[])
	{

		// Find 3 consecutive center-crossing values for I.
		// This determines the value for 360 degrees of phase, and the centerPhase value.
		// Then offset 90 degrees, and see how far off centerPhase the Q value is.
		// Use this as the value of Quadrature Error.
		// CenterPhase takes several iterations to zero in on.

		for (int k=0; k<5; k++)
		{
			// find 1st frequency index where I(index) crosses centerPhase
			for (int i=50; i<900; i++)
			{
				if (Iphase[i] >= centerPhase)
				{
					centerIindex = i;
					break;
				}
			}
			// find 2nd frequency index where I(index) crosses centerPhase
			for (int i=centerIindex+5; i<900; i++)
			{
				if (Iphase[i] <= centerPhase)
				{
					centerI2index = i;
					break;
				}
			}
			
			// find 3rd frequency index where I(index) crosses centerPhase
			for (int i=centerI2index+5; i<1024; i++)
			{
				if (Iphase[i] >= centerPhase)
				{
					centerI3index = i;
					break;
				}
			}

			// Check the results for gross errors that will prevent the algorithm
			// from converging or giving too little range from which to extrapolate results.
			// This can be caused by the test cable being the wrong length.

			if((centerI2index < 50) || (centerI3index < 50))	// The calibration cable is too short
			{
				MessageBox::Show(S"Cannot converge phase calibration.\n\r"
					S"The cable used may be too short.",
					S"Phase Calibration Error", MessageBoxButtons::OK, MessageBoxIcon::Error);
				return(false);
			}

			if(centerI3index < 420)								// The calibration cable is too long
			{
				MessageBox::Show(S"Cannot converge phase calibration.\n\r"
					S"The cable used may be too long.",
					S"Phase Calibration Error", MessageBoxButtons::OK, MessageBoxIcon::Error);
				return(false);
			}

			// Compute level error based on I-I2 and I2-I3 mismatch.
			// Converge on value for centerPhase that equalizes these two distances.

			int error = (centerI3index - centerI2index - (centerI2index - centerIindex)) / 4;
			centerPhase = Iphase[centerI2index+error];
		}

		cxFreqIdx = centerI3index - centerIindex;		// 360 degrees of phase difference as index count
		upperPhase = Iphase[centerI2index-cxFreqIdx/8];	// revise uppercrossover = centerpoint - 45 degrees
		lowerPhase = Iphase[centerI2index+cxFreqIdx/8];	// revise lowercorssover = centerpoint + 45 degrees
		phaseSlope = (float)90.0 / (float)(upperPhase - lowerPhase);	// compute degrees per unit count

		int qerrampl = Qphase[centerIindex + cxFreqIdx/4] - centerPhase;  // Q error as amplitude
		quadratureError = -(float)qerrampl * phaseSlope;						// Q error as degrees

		// Now determine error of measured phase derived from the atan2() function.
		// It's not completely periodic, so generate the error every closest degree.
		// pherror[degrees] = difference.

		for (int k= centerIindex; k<centerI3index; k++)		// sweep over 360 degrees - from +90 to 0 to -180..+180 to -90.
		{
			// derive phase using atan approximation
			double phase = RAD2DEGR * atan2((double)(centerPhase - Qphase[k]),(double)(Iphase[k] - centerPhase));
			// derive phase using cable linearity approximation
			double cphase = 90.0 - 360.0 * (double)(k - centerIindex) / (double)cxFreqIdx;
			if (cphase < -180.0)
				cphase += 360.0;
			double error = cphase - phase;
			if (error < -180.0)
				error += 360.0;
			if (error > 180.0)
				error -= 360.0;
			pherror[(int)cphase+180] = error;	// offset index -->  [-180 .. +180] is the index range
		}

		calibrated = true;
		return(true);
	}
/// Construct amplitude calibration table from raw Magnitude data
bool Detector::AmpCal(int Mag __gc[ , ])
{

	int freq; int level;
	int sumx, sumy, sumxsq, sumysq, sumxy, n;
	int sumf;
	n = 26;											// number of points in regression
    
	//  Use least-squares linear regression to compute  'm' and 'b'  in  "y = mx + b"
	//  'x' is the raw DAC count, 'y' is the value of dbm derived from that raw count. 
	//  Best value for correlation coefficient, 'r' is when the range 0..-25 dbm is used
	//  NB: level in dbm is opposite sign from the array index.
	//  Also compute the level where the curve flatens out (frequency-dependent).