Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: ✨ add support for DRAG #18

Merged
merged 1 commit into from
Jun 30, 2023
Merged

feat: ✨ add support for DRAG #18

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 8 additions & 10 deletions examples/WaveGenDemo/Program.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,27 +42,25 @@ static void Run<T>() where T : unmanaged, IFloatingPointIeee754<T>

var instructions = new List<Instruction>();
var shape = new HannPulseShape();
instructions.Add(new Play(shape, 0, 30e-9, 100e-9, 0.5, 0, 0, ch1));
instructions.Add(new Play(shape, 0, 30e-9, 100e-9, 0.6, 0, 0, ch2));
instructions.Add(new Play(shape, 0, 30e-9, 100e-9, 0.5, 2e-9, 0, 0, ch1));
instructions.Add(new Play(shape, 0, 30e-9, 100e-9, 0.6, 2e-9, 0, 0, ch2));
instructions.Add(new ShiftPhase(0.25 * Math.Tau, ch1));
instructions.Add(new ShiftPhase(-0.25 * Math.Tau, ch2));
instructions.Add(new Play(shape, 200e-9, 30e-9, 100e-9, 0.5, 0, 0, ch1));
instructions.Add(new Play(shape, 200e-9, 30e-9, 100e-9, 0.6, 0, 0, ch2));
instructions.Add(new Play(shape, 200e-9, 30e-9, 100e-9, 0.5, 2e-9, 0, 0, ch1));
instructions.Add(new Play(shape, 200e-9, 30e-9, 100e-9, 0.6, 2e-9, 0, 0, ch2));
instructions.Add(new ShiftFrequency(-100e6, 400e-9, ch1));
instructions.Add(new ShiftFrequency(-250e6, 400e-9, ch2));
instructions.Add(new Play(shape, 400e-9, 30e-9, 100e-9, 0.5, 0, 0, ch1));
instructions.Add(new Play(shape, 400e-9, 30e-9, 100e-9, 0.6, 0, 0, ch2));
instructions.Add(new Play(shape, 400e-9, 200e-9, 0e-9, 0.5, 2e-9, 0, 0, ch1));
instructions.Add(new Play(shape, 400e-9, 200e-9, 0e-9, 0.6, 2e-9, 0, 0, ch2));
instructions.Add(new SetFrequency(0, 600e-9, ch1));
instructions.Add(new SetFrequency(0, 600e-9, ch2));

var tStart = 600e-9;
var count = 0;
while (tStart < 49e-6)
{
//instructions.Add(new Play(shape, tStart, 30e-9, 0, 0.5, 0, 0, ch1));
//instructions.Add(new Play(shape, tStart, 30e-9, 0, 0.6, 0, 0, ch2));
instructions.Add(new Play(shape, tStart, 30e-9, 0e-9, 0.5, 0, 0, ch1));
instructions.Add(new Play(shape, tStart, 30e-9, 0e-9, 0.6, 0, 0, ch2));
instructions.Add(new Play(shape, tStart, 30e-9, 0e-9, 0.5, 2e-9, 0, 0, ch1));
instructions.Add(new Play(shape, tStart, 30e-9, 0e-9, 0.6, 2e-9, 0, 0, ch2));
instructions.Add(new ShiftPhase(0.25 * Math.Tau, ch1));
instructions.Add(new ShiftPhase(-0.25 * Math.Tau, ch2));
tStart += 0.1e-9;
Expand Down
1 change: 1 addition & 0 deletions examples/WaveGenDemo/WaveGenDemo.csproj
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
<OutputType>Exe</OutputType>
<TargetFramework>net7.0</TargetFramework>
<ImplicitUsings>enable</ImplicitUsings>
<ServerGarbageCollection>true</ServerGarbageCollection>
<Nullable>enable</Nullable>
</PropertyGroup>

Expand Down
1 change: 1 addition & 0 deletions src/Qynit.Pulsewave/Play.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ public record Play(
double Width,
double Plateau,
double Amplitude,
double DragCoefficient,
double Frequency,
double Phase,
Channel Channel) : Instruction(nameof(Play), new[] { Channel });
8 changes: 5 additions & 3 deletions src/Qynit.Pulsewave/WaveformGenerator.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -134,10 +134,12 @@ private void Play(Play play)
var width = play.Width;
var plateau = play.Plateau;
var amplitude = play.Amplitude;
var frequency = play.Frequency + context.Frequency + context.FrequencyShift;
var phase = play.Phase + context.Phase + Math.Tau * frequency * tStart;
var dragCoefficient = play.DragCoefficient;
var frameFrequency = context.Frequency + context.FrequencyShift;
var totalFrequency = play.Frequency + frameFrequency;
var phase = (play.Phase + context.Phase + Math.Tau * frameFrequency * tStart) % Math.Tau;
var envelope = new Envelope(pulseShape, width, plateau);
builder.Add(envelope, frequency, tStart, T.CreateChecked(amplitude), T.CreateChecked(phase), T.Zero);
builder.Add(envelope, totalFrequency, tStart, T.CreateChecked(amplitude), T.CreateChecked(phase), T.CreateChecked(dragCoefficient));
}

private class ChannelContext
Expand Down
191 changes: 160 additions & 31 deletions src/Qynit.Pulsewave/WaveformUtils.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
namespace Qynit.Pulsewave;
public static class WaveformUtils
{
public static PooledComplexArray<T> SampleWaveform<T>(EnvelopeInfo envelopeInfo, IPulseShape shape, double width, double plateau)
public static PooledComplexArray<T> SampleEnvelope<T>(EnvelopeInfo envelopeInfo, IPulseShape shape, double width, double plateau)
where T : unmanaged, IFloatingPointIeee754<T>
{
var sampleRate = envelopeInfo.SampleRate;
Expand Down Expand Up @@ -42,6 +42,42 @@ public static PooledComplexArray<T> SampleWaveform<T>(EnvelopeInfo envelopeInfo,
return array;
}

internal static PooledComplexArray<T> SampleWaveform<T>(PulseList<T> pulseList, double sampleRate, int length, int alignLevel) where T : unmanaged, IFloatingPointIeee754<T>
{
var waveform = new PooledComplexArray<T>(length, true);
foreach (var (binKey, bin) in pulseList.Items)
{
foreach (var pulse in bin)
{
var tStart = pulse.Time + pulseList.TimeOffset;
var shape = binKey.Envelope.Shape!;
var width = binKey.Envelope.Width;
var plateau = binKey.Envelope.Plateau;
var frequency = binKey.Frequency;
var iFracStart = TimeAxisUtils.NextFracIndex(tStart, sampleRate, alignLevel);
var iStart = (int)Math.Ceiling(iFracStart);
var envelopeInfo = new EnvelopeInfo(iStart - iFracStart, sampleRate);
using var envelope = SampleEnvelope<T>(envelopeInfo, shape, width, plateau);
var dt = 1 / sampleRate;
var amplitude = pulse.Amplitude * pulseList.AmplitudeMultiplier * IqPair<T>.FromPolarCoordinates(T.One, T.CreateChecked((iStart * dt - tStart) * frequency * Math.Tau));
var complexAmplitude = amplitude.Amplitude;
var dPhase = T.CreateChecked(Math.Tau * frequency * dt);
var arrayIStart = iStart;
var dragAmplitude = amplitude.DragAmplitude;
if (dragAmplitude.I == T.Zero && dragAmplitude.Q == T.Zero)
{
MixAddFrequency(waveform[arrayIStart..], envelope, complexAmplitude, dPhase);
}
else
{
var drag = dragAmplitude * T.CreateChecked(sampleRate);
MixAddFrequencyWithDrag(waveform[arrayIStart..], envelope, complexAmplitude, drag, dPhase);
}
}
}
return waveform;
}

public static void MixAddFrequency<T>(ComplexArraySpan<T> target, ComplexArrayReadOnlySpan<T> source, IqPair<T> amplitude, T dPhase)
where T : unmanaged, IFloatingPointIeee754<T>
{
Expand Down Expand Up @@ -74,56 +110,149 @@ public static void MixAddFrequency<T>(ComplexArraySpan<T> target, ComplexArrayRe
}
var phaserVectorI = new Vector<T>(phaserI);
var phaserVectorQ = new Vector<T>(phaserQ);
var carrierVectorI = phaserVectorI * carrier.I - phaserVectorQ * carrier.Q;
var carrierVectorQ = phaserVectorI * carrier.Q + phaserVectorQ * carrier.I;

for (; i < length - vSize + 1; i += vSize)
{
var sourceVectorI = Unsafe.As<T, Vector<T>>(ref Unsafe.Add(ref sourceI, i));
var sourceVectorQ = Unsafe.As<T, Vector<T>>(ref Unsafe.Add(ref sourceQ, i));
var tempI = sourceVectorI * phaserVectorI - sourceVectorQ * phaserVectorQ;
var tempQ = sourceVectorI * phaserVectorQ + sourceVectorQ * phaserVectorI;
var tempI = sourceVectorI * carrierVectorI - sourceVectorQ * carrierVectorQ;
var tempQ = sourceVectorI * carrierVectorQ + sourceVectorQ * carrierVectorI;
ref var targetVectorI = ref Unsafe.As<T, Vector<T>>(ref Unsafe.Add(ref targetI, i));
ref var targetVectorQ = ref Unsafe.As<T, Vector<T>>(ref Unsafe.Add(ref targetQ, i));
targetVectorI += tempI * carrier.I - tempQ * carrier.Q;
targetVectorQ += tempI * carrier.Q + tempQ * carrier.I;
carrier *= phaserBulk;
targetVectorI += tempI;
targetVectorQ += tempQ;
var newCarrierVectorI = carrierVectorI * phaserBulk.I - carrierVectorQ * phaserBulk.Q;
var newCarrierVectorQ = carrierVectorI * phaserBulk.Q + carrierVectorQ * phaserBulk.I;
carrierVectorI = newCarrierVectorI;
carrierVectorQ = newCarrierVectorQ;
}
carrier = new IqPair<T>(carrierVectorI[0], carrierVectorQ[0]);
}

for (; i < length; i++)
{
var sourceScalarI = Unsafe.Add(ref sourceI, i);
var sourceScalarQ = Unsafe.Add(ref sourceQ, i);
ref var targetScalarI = ref Unsafe.Add(ref targetI, i);
ref var targetScalarQ = ref Unsafe.Add(ref targetQ, i);
targetScalarI += sourceScalarI * carrier.I - sourceScalarQ * carrier.Q;
targetScalarQ += sourceScalarI * carrier.Q + sourceScalarQ * carrier.I;
var sourceIq = new IqPair<T>(Unsafe.Add(ref sourceI, i), Unsafe.Add(ref sourceQ, i));
var targetIq = sourceIq * carrier;
Unsafe.Add(ref targetI, i) += targetIq.I;
Unsafe.Add(ref targetQ, i) += targetIq.Q;
carrier *= phaser;
}
}

internal static PooledComplexArray<T> SampleWaveform<T>(PulseList<T> pulseList, double sampleRate, int length, int alignLevel) where T : unmanaged, IFloatingPointIeee754<T>
public static void MixAddFrequencyWithDrag<T>(ComplexArraySpan<T> target, ComplexArrayReadOnlySpan<T> source, IqPair<T> amplitude, IqPair<T> dragAmplitude, T dPhase)
where T : unmanaged, IFloatingPointIeee754<T>
{
var waveform = new PooledComplexArray<T>(length, true);
foreach (var (binKey, bin) in pulseList.Items)
var length = source.Length;
if (length == 0)
{
foreach (var pulse in bin)
return;
}
Debug.Assert(target.Length >= source.Length);
ref var targetI = ref MemoryMarshal.GetReference(target.DataI);
ref var targetQ = ref MemoryMarshal.GetReference(target.DataQ);
ref var sourceI = ref MemoryMarshal.GetReference(source.DataI);
ref var sourceQ = ref MemoryMarshal.GetReference(source.DataQ);
var carrier = amplitude;
var dragCarrier = dragAmplitude;

// left boundary
{
var sourceIq = new IqPair<T>(Unsafe.Add(ref sourceI, 0), Unsafe.Add(ref sourceQ, 0));
var sourceWithAmplitudeIq = sourceIq * carrier;
if (length == 1)
{
var tStart = pulse.Time + pulseList.TimeOffset;
var shape = binKey.Envelope.Shape!;
var width = binKey.Envelope.Width;
var plateau = binKey.Envelope.Plateau;
var frequency = binKey.Frequency;
var iFracStart = TimeAxisUtils.NextFracIndex(tStart, sampleRate, alignLevel);
var iStart = (int)Math.Ceiling(iFracStart);
var envelopeInfo = new EnvelopeInfo(iStart - iFracStart, sampleRate);
using var envelope = SampleWaveform<T>(envelopeInfo, shape, width, plateau);
var dt = 1 / sampleRate;
var cPhase = pulse.Amplitude.Amplitude * pulseList.AmplitudeMultiplier * IqPair<T>.FromPolarCoordinates(T.One, T.CreateChecked((iStart * dt - tStart) * frequency * Math.Tau));
var dPhase = T.CreateChecked(Math.Tau * frequency * dt);
var arrayIStart = iStart;
MixAddFrequency(waveform[arrayIStart..], envelope, cPhase, dPhase);
Unsafe.Add(ref targetI, 0) += sourceWithAmplitudeIq.I;
Unsafe.Add(ref targetQ, 0) += sourceWithAmplitudeIq.Q;
return;
}
var nextSourceIq = new IqPair<T>(Unsafe.Add(ref sourceI, 1), Unsafe.Add(ref sourceQ, 1));
var diff = nextSourceIq - sourceIq;
var dragWithAmplitudeIq = diff * dragCarrier;
var totalIq = sourceWithAmplitudeIq + dragWithAmplitudeIq;
Unsafe.Add(ref targetI, 0) += totalIq.I;
Unsafe.Add(ref targetQ, 0) += totalIq.Q;
}
var phaser = IqPair<T>.FromPolarCoordinates(T.One, dPhase);
carrier *= phaser;
dragCarrier *= phaser * T.Exp2(-T.One);

var i = 1;
var vSize = Vector<T>.Count;
if (Vector.IsHardwareAccelerated && length >= 2 * vSize + 2)
{
Span<T> phaserI = stackalloc T[vSize];
Span<T> phaserQ = stackalloc T[vSize];
var phaserBulk = IqPair<T>.One;
for (var j = 0; j < vSize; j++)
{
phaserI[j] = phaserBulk.I;
phaserQ[j] = phaserBulk.Q;
phaserBulk *= phaser;
}
var phaserVectorI = new Vector<T>(phaserI);
var phaserVectorQ = new Vector<T>(phaserQ);
var carrierVectorI = phaserVectorI * carrier.I - phaserVectorQ * carrier.Q;
var carrierVectorQ = phaserVectorI * carrier.Q + phaserVectorQ * carrier.I;
var dragCarrierVectorI = phaserVectorI * dragCarrier.I - phaserVectorQ * dragCarrier.Q;
var dragCarrierVectorQ = phaserVectorI * dragCarrier.Q + phaserVectorQ * dragCarrier.I;

for (; i < length - vSize; i += vSize)
{
var sourceVectorI = Unsafe.As<T, Vector<T>>(ref Unsafe.Add(ref sourceI, i));
var sourceVectorQ = Unsafe.As<T, Vector<T>>(ref Unsafe.Add(ref sourceQ, i));
var sourceWithAmplitudeVectorI = sourceVectorI * carrierVectorI - sourceVectorQ * carrierVectorQ;
var sourceWithAmplitudeVectorQ = sourceVectorI * carrierVectorQ + sourceVectorQ * carrierVectorI;
var nextSourceVectorI = Unsafe.As<T, Vector<T>>(ref Unsafe.Add(ref sourceI, i + 1));
var nextSourceVectorQ = Unsafe.As<T, Vector<T>>(ref Unsafe.Add(ref sourceQ, i + 1));
var prevSourceVectorI = Unsafe.As<T, Vector<T>>(ref Unsafe.Add(ref sourceI, i - 1));
var prevSourceVectorQ = Unsafe.As<T, Vector<T>>(ref Unsafe.Add(ref sourceQ, i - 1));
var diffVectorI = nextSourceVectorI - prevSourceVectorI;
var diffVectorQ = nextSourceVectorQ - prevSourceVectorQ;
var dragWithAmplitudeVectorI = diffVectorI * dragCarrierVectorI - diffVectorQ * dragCarrierVectorQ;
var dragWithAmplitudeVectorQ = diffVectorI * dragCarrierVectorQ + diffVectorQ * dragCarrierVectorI;
var totalVectorI = sourceWithAmplitudeVectorI + dragWithAmplitudeVectorI;
var totalVectorQ = sourceWithAmplitudeVectorQ + dragWithAmplitudeVectorQ;
ref var targetVectorI = ref Unsafe.As<T, Vector<T>>(ref Unsafe.Add(ref targetI, i));
ref var targetVectorQ = ref Unsafe.As<T, Vector<T>>(ref Unsafe.Add(ref targetQ, i));
targetVectorI += totalVectorI;
targetVectorQ += totalVectorQ;
var newCarrierVectorI = carrierVectorI * phaserBulk.I - carrierVectorQ * phaserBulk.Q;
var newCarrierVectorQ = carrierVectorI * phaserBulk.Q + carrierVectorQ * phaserBulk.I;
var newDragCarrierVectorI = dragCarrierVectorI * phaserBulk.I - dragCarrierVectorQ * phaserBulk.Q;
var newDragCarrierVectorQ = dragCarrierVectorI * phaserBulk.Q + dragCarrierVectorQ * phaserBulk.I;
carrierVectorI = newCarrierVectorI;
carrierVectorQ = newCarrierVectorQ;
dragCarrierVectorI = newDragCarrierVectorI;
dragCarrierVectorQ = newDragCarrierVectorQ;
}
carrier = new IqPair<T>(carrierVectorI[0], carrierVectorQ[0]);
dragCarrier = new IqPair<T>(dragCarrierVectorI[0], dragCarrierVectorQ[0]);
}

for (; i < length - 1; i++)
{
var sourceIq = new IqPair<T>(Unsafe.Add(ref sourceI, i), Unsafe.Add(ref sourceQ, i));
var nextSourceIq = new IqPair<T>(Unsafe.Add(ref sourceI, i + 1), Unsafe.Add(ref sourceQ, i + 1));
var prevSourceIq = new IqPair<T>(Unsafe.Add(ref sourceI, i - 1), Unsafe.Add(ref sourceQ, i - 1));
var diff = nextSourceIq - prevSourceIq;
var totalIq = sourceIq * carrier + diff * dragCarrier;
Unsafe.Add(ref targetI, i) += totalIq.I;
Unsafe.Add(ref targetQ, i) += totalIq.Q;
carrier *= phaser;
dragCarrier *= phaser;
}

// right boundary
{
dragCarrier *= T.Exp2(T.One);
var sourceIq = new IqPair<T>(Unsafe.Add(ref sourceI, length - 1), Unsafe.Add(ref sourceQ, length - 1));
var prevSourceIq = new IqPair<T>(Unsafe.Add(ref sourceI, length - 2), Unsafe.Add(ref sourceQ, length - 2));
var diff = sourceIq - prevSourceIq;
var totalIq = sourceIq * carrier + diff * dragCarrier;
Unsafe.Add(ref targetI, length - 1) += totalIq.I;
Unsafe.Add(ref targetQ, length - 1) += totalIq.Q;
}
return waveform;
}
}
Loading