The personal website of Scott W Harden
March 28th, 2023

Treemapping with C#

Treemap diagrams display a series of positive numbers using rectangles sized proportional to the value of each number. This page demonstrates how to calculate the size and location of rectangles to create a tree map diagram using C#. Although the following uses System.Drawing to save the tree map as a Bitmap image, these concepts may be combined with information on the C# Data Visualization page to create treemap diagrams using SkiaSharp, WPF, or other graphics technologies.

The tree map above was generated from random data using the following C# code:

// Create sample data. Data must be sorted large to small.
double[] sortedValues = Enumerable.Range(0, 40)
    .Select(x => (double)Random.Shared.Next(10, 100))
    .OrderByDescending(x => x)
    .ToArray();

// Create an array of labels in the same order as the sorted data.
string[] labels = sortedValues.Select(x => x.ToString()).ToArray();

// Calculate the size and position of all rectangles in the tree map
int width = 600;
int height = 400;
RectangleF[] rectangles = TreeMap.GetRectangles(sortedValues, width, height);

// Create an image to draw on (with 1px extra to make room for the outline)
using Bitmap bmp = new(width + 1, height + 1);
using Graphics gfx = Graphics.FromImage(bmp);
using Font fnt = new("Consolas", 8);
using SolidBrush brush = new(Color.Black);
gfx.Clear(Color.White);

// Draw and label each rectangle
for (int i = 0; i < rectangles.Length; i++)
{
    brush.Color = Color.FromArgb(
        red: Random.Shared.Next(150, 250),
        green: Random.Shared.Next(150, 250),
        blue: Random.Shared.Next(150, 250));

    gfx.FillRectangle(brush, rectangles[i]);
    gfx.DrawRectangle(Pens.Black, rectangles[i]);
    gfx.DrawString(labels[i], fnt, Brushes.Black, rectangles[i].X, rectangles[i].Y);
}

// Save the output
bmp.Save("treemap.bmp");

Treemap Logic

The previous code block focuses on data generation and display, but hides the tree map calculations behind the TreeMap class. Below is the code for that class. It is self-contained static class and exposes a single static method which takes a pre-sorted array of values and returns tree map rectangles ready to display on an image.

💡 Although the System.Drawing.Common is a Windows-only library (as of .NET 7), System.Drawing.Primitives is a cross-platform package that provides the RectangleF structure used in the tree map class. See the SkiaSharp Quickstart to learn how to create image files using cross-platform .NET code.

public static class TreeMap
{
    public static RectangleF[] GetRectangles(double[] values, int width, int height)
    {
        for (int i = 1; i < values.Length; i++)
            if (values[i] > values[i - 1])
                throw new ArgumentException("values must be ordered large to small");

        var slice = GetSlice(values, 1, 0.35);
        var rectangles = GetRectangles(slice, width, height);
        return rectangles.Select(x => x.ToRectF()).ToArray();
    }

    private class Slice
    {
        public double Size { get; }
        public IEnumerable<double> Values { get; }
        public Slice[] Children { get; }

        public Slice(double size, IEnumerable<double> values, Slice sub1, Slice sub2)
        {
            Size = size;
            Values = values;
            Children = new Slice[] { sub1, sub2 };
        }

        public Slice(double size, double finalValue)
        {
            Size = size;
            Values = new double[] { finalValue };
            Children = Array.Empty<Slice>();
        }
    }

    private class SliceResult
    {
        public double ElementsSize { get; }
        public IEnumerable<double> Elements { get; }
        public IEnumerable<double> RemainingElements { get; }

        public SliceResult(double elementsSize, IEnumerable<double> elements, IEnumerable<double> remainingElements)
        {
            ElementsSize = elementsSize;
            Elements = elements;
            RemainingElements = remainingElements;
        }
    }

    private class SliceRectangle
    {
        public Slice Slice { get; set; }
        public float X { get; set; }
        public float Y { get; set; }
        public float Width { get; set; }
        public float Height { get; set; }
        public SliceRectangle(Slice slice) => Slice = slice;
        public RectangleF ToRectF() => new(X, Y, Width, Height);
    }

    private static Slice GetSlice(IEnumerable<double> elements, double totalSize, double sliceWidth)
    {
        if (elements.Count() == 1)
            return new Slice(totalSize, elements.Single());

        SliceResult sr = GetElementsForSlice(elements, sliceWidth);
        Slice child1 = GetSlice(sr.Elements, sr.ElementsSize, sliceWidth);
        Slice child2 = GetSlice(sr.RemainingElements, 1 - sr.ElementsSize, sliceWidth);
        return new Slice(totalSize, elements, child1, child2);
    }

    private static SliceResult GetElementsForSlice(IEnumerable<double> elements, double sliceWidth)
    {
        var elementsInSlice = new List<double>();
        var remainingElements = new List<double>();
        double current = 0;
        double total = elements.Sum();

        foreach (var element in elements)
        {
            if (current > sliceWidth)
                remainingElements.Add(element);
            else
            {
                elementsInSlice.Add(element);
                current += element / total;
            }
        }

        return new SliceResult(current, elementsInSlice, remainingElements);
    }

    private static IEnumerable<SliceRectangle> GetRectangles(Slice slice, int width, int height)
    {
        SliceRectangle area = new(slice) { Width = width, Height = height };

        foreach (var rect in GetRectangles(area))
        {
            if (rect.X + rect.Width > area.Width)
                rect.Width = area.Width - rect.X;

            if (rect.Y + rect.Height > area.Height)
                rect.Height = area.Height - rect.Y;

            yield return rect;
        }
    }

    private static IEnumerable<SliceRectangle> GetRectangles(SliceRectangle sliceRectangle)
    {
        var isHorizontalSplit = sliceRectangle.Width >= sliceRectangle.Height;
        var currentPos = 0;
        foreach (var subSlice in sliceRectangle.Slice.Children)
        {
            var subRect = new SliceRectangle(subSlice);
            int rectSize;

            if (isHorizontalSplit)
            {
                rectSize = (int)Math.Round(sliceRectangle.Width * subSlice.Size);
                subRect.X = sliceRectangle.X + currentPos;
                subRect.Y = sliceRectangle.Y;
                subRect.Width = rectSize;
                subRect.Height = sliceRectangle.Height;
            }
            else
            {
                rectSize = (int)Math.Round(sliceRectangle.Height * subSlice.Size);
                subRect.X = sliceRectangle.X;
                subRect.Y = sliceRectangle.Y + currentPos;
                subRect.Width = sliceRectangle.Width;
                subRect.Height = rectSize;
            }

            currentPos += rectSize;

            if (subSlice.Values.Count() > 1)
            {
                foreach (var sr in GetRectangles(subRect))
                {
                    yield return sr;
                }
            }
            else if (subSlice.Values.Count() == 1)
            {
                yield return subRect;
            }
        }
    }
}

Source Code Complexity Analysis

A few days ago I wrote an article describing how to programmatically generate .NET source code analytics using C#. Using these tools I analyzed the source code for all classes in a large project (ScottPlot.NET). The following tree map displays every class in the project as a rectangle sized according to number of lines of code and colored according to maintainability.

In this diagram large rectangles represent classes with the most code, and red color indicates classes that are difficult to maintain.

💡 I'm using a perceptually uniform colormap (similar to Turbo)provided by the ScottPlot provided by the ScottPlot NuGet package. See ScottPlot's colormaps gallery for all available colormaps.

💡 The Maintainability Index is a value between 0 (worst) and 100 (best) that represents the relative ease of maintaining the code. It's calculated from a combination of Halstead complexity (size of the compiled code), Cyclomatic complexity (number of paths that can be taken through the code), and the total number of lines of code.

Conclusions

  • Generation of tree map diagrams can be achieved using recursive programming

  • The static class above makes it easy to generate tree maps in C#

  • ScottPlot's AxisTicksRender class may be difficult to maintain

References

Markdown source code last modified on March 8th, 2023
---
Title: Treemapping with C#
Description: How to create a treemap diagram using C#
Date: 2023-03-28 00:32AM EST
Tags: csharp, graphics
---

# Treemapping with C# 

**Treemap diagrams display a series of positive numbers using rectangles sized proportional to the value of each number.** This page demonstrates how to calculate the size and location of rectangles to create a tree map diagram using C#. Although the following uses System.Drawing to save the tree map as a Bitmap image, these concepts may be combined with information on the [C# Data Visualization](https://swharden.com/csdv/) page to create treemap diagrams using SkiaSharp, WPF, or other graphics technologies.

<img src="treemap.png" class="img-fluid d-block mx-auto shadow my-5">

The tree map above was generated from random data using the following C# code:

```cs
// Create sample data. Data must be sorted large to small.
double[] sortedValues = Enumerable.Range(0, 40)
    .Select(x => (double)Random.Shared.Next(10, 100))
    .OrderByDescending(x => x)
    .ToArray();

// Create an array of labels in the same order as the sorted data.
string[] labels = sortedValues.Select(x => x.ToString()).ToArray();

// Calculate the size and position of all rectangles in the tree map
int width = 600;
int height = 400;
RectangleF[] rectangles = TreeMap.GetRectangles(sortedValues, width, height);

// Create an image to draw on (with 1px extra to make room for the outline)
using Bitmap bmp = new(width + 1, height + 1);
using Graphics gfx = Graphics.FromImage(bmp);
using Font fnt = new("Consolas", 8);
using SolidBrush brush = new(Color.Black);
gfx.Clear(Color.White);

// Draw and label each rectangle
for (int i = 0; i < rectangles.Length; i++)
{
    brush.Color = Color.FromArgb(
        red: Random.Shared.Next(150, 250),
        green: Random.Shared.Next(150, 250),
        blue: Random.Shared.Next(150, 250));

    gfx.FillRectangle(brush, rectangles[i]);
    gfx.DrawRectangle(Pens.Black, rectangles[i]);
    gfx.DrawString(labels[i], fnt, Brushes.Black, rectangles[i].X, rectangles[i].Y);
}

// Save the output
bmp.Save("treemap.bmp");
```

## Treemap Logic

The previous code block focuses on data generation and display, but hides the tree map calculations behind the `TreeMap` class. Below is the code for that class. It is self-contained static class and exposes a single static method which takes a pre-sorted array of values and returns tree map rectangles ready to display on an image.

> 💡 Although the `System.Drawing.Common` is a Windows-only library ([as of .NET 7](https://github.com/dotnet/designs/blob/main/accepted/2021/system-drawing-win-only/system-drawing-win-only.md)), `System.Drawing.Primitives` is a cross-platform package that provides the `RectangleF` structure used in the tree map class. See the [SkiaSharp Quickstart](https://swharden.com/csdv/skiasharp/quickstart-console/) to learn how to create image files using cross-platform .NET code.

```cs
public static class TreeMap
{
    public static RectangleF[] GetRectangles(double[] values, int width, int height)
    {
        for (int i = 1; i < values.Length; i++)
            if (values[i] > values[i - 1])
                throw new ArgumentException("values must be ordered large to small");

        var slice = GetSlice(values, 1, 0.35);
        var rectangles = GetRectangles(slice, width, height);
        return rectangles.Select(x => x.ToRectF()).ToArray();
    }

    private class Slice
    {
        public double Size { get; }
        public IEnumerable<double> Values { get; }
        public Slice[] Children { get; }

        public Slice(double size, IEnumerable<double> values, Slice sub1, Slice sub2)
        {
            Size = size;
            Values = values;
            Children = new Slice[] { sub1, sub2 };
        }

        public Slice(double size, double finalValue)
        {
            Size = size;
            Values = new double[] { finalValue };
            Children = Array.Empty<Slice>();
        }
    }

    private class SliceResult
    {
        public double ElementsSize { get; }
        public IEnumerable<double> Elements { get; }
        public IEnumerable<double> RemainingElements { get; }

        public SliceResult(double elementsSize, IEnumerable<double> elements, IEnumerable<double> remainingElements)
        {
            ElementsSize = elementsSize;
            Elements = elements;
            RemainingElements = remainingElements;
        }
    }

    private class SliceRectangle
    {
        public Slice Slice { get; set; }
        public float X { get; set; }
        public float Y { get; set; }
        public float Width { get; set; }
        public float Height { get; set; }
        public SliceRectangle(Slice slice) => Slice = slice;
        public RectangleF ToRectF() => new(X, Y, Width, Height);
    }

    private static Slice GetSlice(IEnumerable<double> elements, double totalSize, double sliceWidth)
    {
        if (elements.Count() == 1)
            return new Slice(totalSize, elements.Single());

        SliceResult sr = GetElementsForSlice(elements, sliceWidth);
        Slice child1 = GetSlice(sr.Elements, sr.ElementsSize, sliceWidth);
        Slice child2 = GetSlice(sr.RemainingElements, 1 - sr.ElementsSize, sliceWidth);
        return new Slice(totalSize, elements, child1, child2);
    }

    private static SliceResult GetElementsForSlice(IEnumerable<double> elements, double sliceWidth)
    {
        var elementsInSlice = new List<double>();
        var remainingElements = new List<double>();
        double current = 0;
        double total = elements.Sum();

        foreach (var element in elements)
        {
            if (current > sliceWidth)
                remainingElements.Add(element);
            else
            {
                elementsInSlice.Add(element);
                current += element / total;
            }
        }

        return new SliceResult(current, elementsInSlice, remainingElements);
    }

    private static IEnumerable<SliceRectangle> GetRectangles(Slice slice, int width, int height)
    {
        SliceRectangle area = new(slice) { Width = width, Height = height };

        foreach (var rect in GetRectangles(area))
        {
            if (rect.X + rect.Width > area.Width)
                rect.Width = area.Width - rect.X;

            if (rect.Y + rect.Height > area.Height)
                rect.Height = area.Height - rect.Y;

            yield return rect;
        }
    }

    private static IEnumerable<SliceRectangle> GetRectangles(SliceRectangle sliceRectangle)
    {
        var isHorizontalSplit = sliceRectangle.Width >= sliceRectangle.Height;
        var currentPos = 0;
        foreach (var subSlice in sliceRectangle.Slice.Children)
        {
            var subRect = new SliceRectangle(subSlice);
            int rectSize;

            if (isHorizontalSplit)
            {
                rectSize = (int)Math.Round(sliceRectangle.Width * subSlice.Size);
                subRect.X = sliceRectangle.X + currentPos;
                subRect.Y = sliceRectangle.Y;
                subRect.Width = rectSize;
                subRect.Height = sliceRectangle.Height;
            }
            else
            {
                rectSize = (int)Math.Round(sliceRectangle.Height * subSlice.Size);
                subRect.X = sliceRectangle.X;
                subRect.Y = sliceRectangle.Y + currentPos;
                subRect.Width = sliceRectangle.Width;
                subRect.Height = rectSize;
            }

            currentPos += rectSize;

            if (subSlice.Values.Count() > 1)
            {
                foreach (var sr in GetRectangles(subRect))
                {
                    yield return sr;
                }
            }
            else if (subSlice.Values.Count() == 1)
            {
                yield return subRect;
            }
        }
    }
}
```

## Source Code Complexity Analysis

A few days ago I wrote an article describing how to [programmatically generate .NET source code analytics using C#](https://swharden.com/blog/2023-03-05-dotnet-code-analysis/). Using these tools I analyzed the source code for all classes in a large project ([ScottPlot.NET](https://scottplot.net)). The following tree map displays every class in the project as a rectangle sized according to number of lines of code and colored according to maintainability. 

<a href="code-report.png"><img src="code-report.png" class="img-fluid d-inline-block mx-auto shadow mt-5"></a>

<img src="turbo.png" class="img-fluid d-inline-block mx-auto shadow mb-5">

In this diagram large rectangles represent classes with the most code, and red color indicates classes that are difficult to maintain. 

> 💡 I'm using a perceptually uniform colormap (similar to [Turbo](https://ai.googleblog.com/2019/08/turbo-improved-rainbow-colormap-for.html))provided by the [ScottPlot](https://scottplot.net) provided by the [ScottPlot](https://scottplot.net) NuGet package. See ScottPlot's [colormaps gallery](https://scottplot.net/cookbook/4.1/colormaps/) for all available colormaps.

> 💡 The [Maintainability Index](https://learn.microsoft.com/en-us/visualstudio/code-quality/code-metrics-maintainability-index-range-and-meaning) is a value between 0 (worst) and 100 (best) that represents the relative ease of maintaining the code. It's calculated from a combination of [Halstead complexity](https://en.wikipedia.org/wiki/Halstead_complexity_measures) (size of the compiled code), [Cyclomatic complexity](https://en.wikipedia.org/wiki/Cyclomatic_complexity) (number of paths that can be taken through the code), and the total number of lines of code.

## Conclusions

* Generation of tree map diagrams can be achieved using recursive programming

* The static class above makes it easy to generate tree maps in C#

* ScottPlot's AxisTicksRender class may be difficult to maintain

## References

* This blog post is spillover from ScottPlot issues [#1479](https://github.com/ScottPlot/ScottPlot/issues/1479) and [#2454](https://github.com/ScottPlot/ScottPlot/issues/2454).

* Code here was heavily influenced by [The Never Ending Journey](http://pascallaurin42.blogspot.com/2013/12/implementing-treemap-in-c.html) (Dec 29, 2013)

* [Treemapping](https://en.wikipedia.org/wiki/Treemapping) (Wikipedia)

* [D3 Treemap](https://d3-graph-gallery.com/treemap.html)

* [Squarified Treemaps](https://www.win.tue.nl/~vanwijk/stm.pdf) (Bruls et al.)

* StackOverflow question [32548949](https://stackoverflow.com/questions/32548949/from-c-sharp-serverside-is-there-anyway-to-generate-a-treemap-and-save-as-an-im/37154938#37154938)

* [C# Data Visualization](https://swharden.com/csdv/)
March 5th, 2023

Generate .NET Code Metrics from Console Applications

This page describes how to use the Microsoft.CodeAnalysis.Metrics package to perform source code analysis of .NET assemblies from a console application. Visual Studio users can perform source code analysis by clicking the "Analyze" dropdown menu and selecting "Calculate Code Metrics", but I sought to automate this process so I can generate custom code analysis reports from console applications as part of my CI pipeline.

Performing Code Analysis

Step 1: Add the Microsoft.CodeAnalysis.Metrics package to your project:

dotnet add package Microsoft.CodeAnalysis.Metrics

Step 2: Perform code analysis:

dotnet build -target:Metrics

Note that multi-targeted projects must append --framework net6.0 to specify a single platform target to use for code analysis.

Step 3: Inspect analysis results in ProjectName.Metrics.xml

<?xml version="1.0" encoding="utf-8"?>
<CodeMetricsReport Version="1.0">
  <Targets>
    <Target Name="ScottPlot.csproj">
      <Assembly Name="ScottPlot, Version=4.1.61.0">
        <Metrics>
          <Metric Name="MaintainabilityIndex" Value="81" />
          <Metric Name="CyclomaticComplexity" Value="6324" />
          <Metric Name="ClassCoupling" Value="664" />
          <Metric Name="DepthOfInheritance" Value="3" />
          <Metric Name="SourceLines" Value="35360" />
          <Metric Name="ExecutableLines" Value="10208" />
        </Metrics>
        <Namespaces>
        ...

Parsing the Analysis XML File

The code analysis XML contains information about every assembly, namespace, type, and function in the whole code base! There is a lot of possible information to extract, but the code below is enough to get us started extracting basic metric information for every type in the code base.

using System;
using System.IO;
using System.Linq;
using System.Text;
using System.Xml.Linq;
using System.Collections.Generic;

/// <summary>
/// Display a particular metric for every type in an assembly.
/// </summary>
void RankTypes(string xmlFilePath, string metricName = "CyclomaticComplexity", bool highToLow = true)
{
    string xmlText = File.ReadAllText(xmlFilePath);
    XDocument doc = XDocument.Parse(xmlText);
    XElement assembly = doc.Descendants("Assembly").First();

    var rankedTypes = GetMetricByType(assembly, metricName).OrderBy(x => x.Value).ToArray();
    if (highToLow)
        Array.Reverse(rankedTypes);

    Console.WriteLine($"Types ranked by {metricName}:");
    foreach (var type in rankedTypes)
        Console.WriteLine($"{type.Value:N0}\t{type.Key}");
}

Dictionary<string, int> GetMetricByType(XElement assembly, string metricName)
{
    Dictionary<string, int> metricByType = new();

    foreach (XElement namespaceElement in assembly.Element("Namespaces")!.Elements("Namespace"))
    {
        foreach (XElement namedType in namespaceElement.Elements("Types").Elements("NamedType"))
        {
            XElement metric = namedType.Element("Metrics")!.Elements("Metric")
                .Where(x => x.Attribute("Name")!.Value == metricName)
                .Single();
            string typeName = namedType.Attribute("Name")!.Value;
            string namespaceName = namespaceElement.Attribute("Name")!.Value;
            string fullTypeName = $"{namespaceName}.{typeName}";
            metricByType[fullTypeName] = int.Parse(metric.Attribute("Value")!.Value!.ToString());
        }
    }

    return metricByType;
}

Querying Code Analysis Results

Specific metrics of interest will vary, but here are some code examples demonstrating how to parse the code metrics file to display useful information. For these examples I run the code analysis command above to generate ScottPlot.Metrics.xml from the ScottPlot code base and use the code above to generate various reports.

Rank Types by Cyclomatic Complexity

Cyclomatic complexity is a measure of the number of different paths that can be taken through a computer program, and it is often used as an indicator for difficult-to-maintain code. Some CI systems even prevent the merging of pull requests if their cyclomatic complexity exceeds a predefined threshold! Although I don't intend to gate pull requests by complexity at this time, I would like to gain insight into which classes are the most complex as a way to quantitatively target my code maintenance and efforts.

RankTypes("ScottPlot.Metrics.xml", "CyclomaticComplexity");
517     ScottPlot.Plot
218     ScottPlot.Plottable.SignalPlotBase<T>
173     ScottPlot.Plottable.ScatterPlot
139     ScottPlot.Settings
120     ScottPlot.Ticks.TickCollection
118     ScottPlot.Renderable.Axis
114     ScottPlot.Drawing.Colormap
113     ScottPlot.Control.ControlBackEnd
109     ScottPlot.DataGen
99      ScottPlot.Plottable.AxisLineVector
98      ScottPlot.Plottable.Heatmap
95      ScottPlot.Tools
93      ScottPlot.Plottable.RepeatingAxisLine
91      ScottPlot.Plottable.PopulationPlot
85      ScottPlot.Plottable.AxisLine
83      ScottPlot.Plottable.AxisSpan
77      ScottPlot.Plottable.RadialGaugePlot
...

Rank Types by Lines of Code

Similarly, ranking all my project's types by how many lines of code they contain can give me insight into which types may benefit most from refactoring.

RankTypes("ScottPlot.Metrics.xml", "SourceLines");
Types ranked by SourceLines:
4,155   ScottPlot.Plot
1,182   ScottPlot.DataGen
954     ScottPlot.Plottable.SignalPlotBase<T>
726     ScottPlot.Control.ControlBackEnd
670     ScottPlot.Ticks.TickCollection
670     ScottPlot.Settings
630     ScottPlot.Plottable.ScatterPlot
600     ScottPlot.Renderable.Axis
477     ScottPlot.Statistics.Common
454     ScottPlot.Tools
451     ScottPlot.Plottable.PopulationPlot
432     ScottPlot.Drawing.GDI
343     ScottPlot.Plottable.SignalPlotXYGeneric<TX, TY>
336     ScottPlot.Plottable.RepeatingAxisLine
335     ScottPlot.Drawing.Colormap
332     ScottPlot.Plottable.AxisLineVector
...

Rank Types by Maintainability

The Maintainability Index is a value between 0 (worst) and 100 (best) that represents the relative ease of maintaining the code. It's calculated from a combination of Halstead complexity (size of the compiled code), Cyclomatic complexity (number of paths that can be taken through the code), and the total number of lines of code.

MaintainabilityIndex = 171 
  - 5.2 * Math.Log(HalsteadVolume) 
  - 0.23 * CyclomaticComplexity
  - 16.2 * Math.Log(LinesOfCCode);

The maintainability index is calculated by Microsoft.CodeAnalysis.Metrics so we don't have to. I don't know how Microsoft arrived at their weights for this formula, but the overall idea is described here.

RankTypes("ScottPlot.Metrics.xml", "MaintainabilityIndex", highToLow: false);
43      ScottPlot.Drawing.Tools
48      ScottPlot.Statistics.Interpolation.Cubic
48      ScottPlot.Statistics.Interpolation.PeriodicSpline
49      ScottPlot.Statistics.Interpolation.EndSlopeSpline
49      ScottPlot.Statistics.Interpolation.NaturalSpline
50      ScottPlot.Renderable.AxisTicksRender
54      ScottPlot.Statistics.Interpolation.CatmullRom
55      ScottPlot.Statistics.Interpolation.SplineInterpolator
57      ScottPlot.DataGen
58      ScottPlot.DataStructures.SegmentedTree<T>
58      ScottPlot.MarkerShapes.Hashtag
58      ScottPlot.Ticks.TickCollection
59      ScottPlot.MarkerShapes.Asterisk
59      ScottPlot.Plottable.SignalPlotXYGeneric<TX, TY>
59      ScottPlot.Statistics.Interpolation.Bezier
60      ScottPlot.Statistics.Interpolation.Chaikin
61      ScottPlot.Generate
61      ScottPlot.Plot
61      ScottPlot.Statistics.Finance
...

Create Custom HTML Reports

With a little more effort you can generate HTML reports that use tables and headings to highlight useful code metrics and draw attention to types that could benefit from refactoring to improve maintainability.

Conclusions

Microsoft's official Microsoft.CodeAnalysis.Metrics NuGet package is a useful tool for analyzing assemblies, navigating through namespaces, types, properties, and methods, and evaluating their metrics. Since these analyses can be performed using console applications, they can be easily integrated into CI pipelines or used to create standalone code analysis applications. Future projects can build on the concepts described here to create graphical visualizations of code metrics in large projects.

Resources

Markdown source code last modified on March 5th, 2023
---
Title: .NET Source Code Analysis
Description: How to analyze source code metrics of .NET assemblies from a console application
Date: 2023-03-05 3:20PM EST
Tags: csharp
---

# Generate .NET Code Metrics from Console Applications

**This page describes how to use the [Microsoft.CodeAnalysis.Metrics](https://www.nuget.org/packages/Microsoft.CodeAnalysis.Metrics/) package to perform source code analysis of .NET assemblies from a console application.** Visual Studio users can perform source code analysis by clicking the "Analyze" dropdown menu and selecting "Calculate Code Metrics", but I sought to automate this process so I can generate custom code analysis reports from console applications as part of my CI pipeline.

## Performing Code Analysis

**Step 1:** Add the [`Microsoft.CodeAnalysis.Metrics`](https://www.nuget.org/packages/Microsoft.CodeAnalysis.Metrics/) package to your project:

```bash
dotnet add package Microsoft.CodeAnalysis.Metrics
```

**Step 2:** Perform code analysis:

```bash
dotnet build -target:Metrics
```

Note that multi-targeted projects must append `--framework net6.0` to specify a single platform target to use for code analysis.

**Step 3:** Inspect analysis results in `ProjectName.Metrics.xml`

```xml
<?xml version="1.0" encoding="utf-8"?>
<CodeMetricsReport Version="1.0">
  <Targets>
    <Target Name="ScottPlot.csproj">
      <Assembly Name="ScottPlot, Version=4.1.61.0">
        <Metrics>
          <Metric Name="MaintainabilityIndex" Value="81" />
          <Metric Name="CyclomaticComplexity" Value="6324" />
          <Metric Name="ClassCoupling" Value="664" />
          <Metric Name="DepthOfInheritance" Value="3" />
          <Metric Name="SourceLines" Value="35360" />
          <Metric Name="ExecutableLines" Value="10208" />
        </Metrics>
        <Namespaces>
        ...
```

## Parsing the Analysis XML File

The code analysis XML contains information about every assembly, namespace, type, and function in the whole code base! There is a lot of possible information to extract, but the code below is enough to get us started extracting basic metric information for every type in the code base.

```cs
using System;
using System.IO;
using System.Linq;
using System.Text;
using System.Xml.Linq;
using System.Collections.Generic;

/// <summary>
/// Display a particular metric for every type in an assembly.
/// </summary>
void RankTypes(string xmlFilePath, string metricName = "CyclomaticComplexity", bool highToLow = true)
{
    string xmlText = File.ReadAllText(xmlFilePath);
    XDocument doc = XDocument.Parse(xmlText);
    XElement assembly = doc.Descendants("Assembly").First();

    var rankedTypes = GetMetricByType(assembly, metricName).OrderBy(x => x.Value).ToArray();
    if (highToLow)
        Array.Reverse(rankedTypes);

    Console.WriteLine($"Types ranked by {metricName}:");
    foreach (var type in rankedTypes)
        Console.WriteLine($"{type.Value:N0}\t{type.Key}");
}

Dictionary<string, int> GetMetricByType(XElement assembly, string metricName)
{
    Dictionary<string, int> metricByType = new();

    foreach (XElement namespaceElement in assembly.Element("Namespaces")!.Elements("Namespace"))
    {
        foreach (XElement namedType in namespaceElement.Elements("Types").Elements("NamedType"))
        {
            XElement metric = namedType.Element("Metrics")!.Elements("Metric")
                .Where(x => x.Attribute("Name")!.Value == metricName)
                .Single();
            string typeName = namedType.Attribute("Name")!.Value;
            string namespaceName = namespaceElement.Attribute("Name")!.Value;
            string fullTypeName = $"{namespaceName}.{typeName}";
            metricByType[fullTypeName] = int.Parse(metric.Attribute("Value")!.Value!.ToString());
        }
    }

    return metricByType;
}
```

## Querying Code Analysis Results

**Specific metrics of interest will vary, but here are some code examples demonstrating how to parse the code metrics file to display useful information.** For these examples I run the code analysis command above to generate [ScottPlot.Metrics.xml](ScottPlot.Metrics.xml.zip) from the [ScottPlot](https://scottplot.net) code base and use the code above to generate various reports.

### Rank Types by Cyclomatic Complexity

[Cyclomatic complexity](https://en.wikipedia.org/wiki/Cyclomatic_complexity) is a measure of the number of different paths that can be taken through a computer program, and it is often used as an indicator for difficult-to-maintain code. Some CI systems even prevent the merging of pull requests if their cyclomatic complexity exceeds a predefined threshold! Although I don't intend to gate pull requests by complexity at this time, I would like to gain insight into which classes are the most complex as a way to quantitatively target my code maintenance and efforts.

```cs
RankTypes("ScottPlot.Metrics.xml", "CyclomaticComplexity");
```

```txt
517     ScottPlot.Plot
218     ScottPlot.Plottable.SignalPlotBase<T>
173     ScottPlot.Plottable.ScatterPlot
139     ScottPlot.Settings
120     ScottPlot.Ticks.TickCollection
118     ScottPlot.Renderable.Axis
114     ScottPlot.Drawing.Colormap
113     ScottPlot.Control.ControlBackEnd
109     ScottPlot.DataGen
99      ScottPlot.Plottable.AxisLineVector
98      ScottPlot.Plottable.Heatmap
95      ScottPlot.Tools
93      ScottPlot.Plottable.RepeatingAxisLine
91      ScottPlot.Plottable.PopulationPlot
85      ScottPlot.Plottable.AxisLine
83      ScottPlot.Plottable.AxisSpan
77      ScottPlot.Plottable.RadialGaugePlot
...
```

### Rank Types by Lines of Code
Similarly, ranking all my project's types by how many lines of code they contain can give me insight into which types may benefit most from refactoring.

```cs
RankTypes("ScottPlot.Metrics.xml", "SourceLines");
```

```txt
Types ranked by SourceLines:
4,155   ScottPlot.Plot
1,182   ScottPlot.DataGen
954     ScottPlot.Plottable.SignalPlotBase<T>
726     ScottPlot.Control.ControlBackEnd
670     ScottPlot.Ticks.TickCollection
670     ScottPlot.Settings
630     ScottPlot.Plottable.ScatterPlot
600     ScottPlot.Renderable.Axis
477     ScottPlot.Statistics.Common
454     ScottPlot.Tools
451     ScottPlot.Plottable.PopulationPlot
432     ScottPlot.Drawing.GDI
343     ScottPlot.Plottable.SignalPlotXYGeneric<TX, TY>
336     ScottPlot.Plottable.RepeatingAxisLine
335     ScottPlot.Drawing.Colormap
332     ScottPlot.Plottable.AxisLineVector
...
```

### Rank Types by Maintainability

The [Maintainability Index](https://learn.microsoft.com/en-us/visualstudio/code-quality/code-metrics-maintainability-index-range-and-meaning) is a value between 0 (worst) and 100 (best) that represents the relative ease of maintaining the code. It's calculated from a combination of [Halstead complexity](https://en.wikipedia.org/wiki/Halstead_complexity_measures) (size of the compiled code), [Cyclomatic complexity](https://en.wikipedia.org/wiki/Cyclomatic_complexity) (number of paths that can be taken through the code), and the total number of lines of code.

```cs
MaintainabilityIndex = 171 
  - 5.2 * Math.Log(HalsteadVolume) 
  - 0.23 * CyclomaticComplexity
  - 16.2 * Math.Log(LinesOfCCode);
```

The maintainability index is calculated by `Microsoft.CodeAnalysis.Metrics` so we don't have to. I don't know how Microsoft arrived at their weights for this formula, but the overall idea is described [here](https://learn.microsoft.com/en-us/visualstudio/code-quality/code-metrics-maintainability-index-range-and-meaning).

```cs
RankTypes("ScottPlot.Metrics.xml", "MaintainabilityIndex", highToLow: false);
```

```txt
43      ScottPlot.Drawing.Tools
48      ScottPlot.Statistics.Interpolation.Cubic
48      ScottPlot.Statistics.Interpolation.PeriodicSpline
49      ScottPlot.Statistics.Interpolation.EndSlopeSpline
49      ScottPlot.Statistics.Interpolation.NaturalSpline
50      ScottPlot.Renderable.AxisTicksRender
54      ScottPlot.Statistics.Interpolation.CatmullRom
55      ScottPlot.Statistics.Interpolation.SplineInterpolator
57      ScottPlot.DataGen
58      ScottPlot.DataStructures.SegmentedTree<T>
58      ScottPlot.MarkerShapes.Hashtag
58      ScottPlot.Ticks.TickCollection
59      ScottPlot.MarkerShapes.Asterisk
59      ScottPlot.Plottable.SignalPlotXYGeneric<TX, TY>
59      ScottPlot.Statistics.Interpolation.Bezier
60      ScottPlot.Statistics.Interpolation.Chaikin
61      ScottPlot.Generate
61      ScottPlot.Plot
61      ScottPlot.Statistics.Finance
...
```

### Create Custom HTML Reports

With a little more effort you can generate HTML reports that use tables and headings to highlight useful code metrics and draw attention to types that could benefit from refactoring to improve maintainability.

* View the sample report: [report.html](report.html)
* Download the code used to generate it: [CodeAnalysisReport.zip](CodeAnalysisReport.zip)

<a href="report.html"><img src="report.png" class="img-fluid d-block my-5 border shadow mx-auto"></a>

## Conclusions

Microsoft's official [Microsoft.CodeAnalysis.Metrics](https://www.nuget.org/packages/Microsoft.CodeAnalysis.Metrics/) NuGet package is a useful tool for analyzing assemblies, navigating through namespaces, types, properties, and methods, and evaluating their metrics. Since these analyses can be performed using console applications, they can be easily integrated into CI pipelines or used to create standalone code analysis applications. Future projects can build on the concepts described here to create graphical visualizations of code metrics in large projects.

## Resources

* [Code metrics values](https://learn.microsoft.com/en-us/visualstudio/code-quality/code-metrics-values) - official documentation of the code metrics analysis system

* [Visual Studio source code analysis](https://learn.microsoft.com/en-us/visualstudio/code-quality/roslyn-analyzers-overview)

* [Microsoft.CodeAnalysis.Metrics NuGet Package](https://www.nuget.org/packages/Microsoft.CodeAnalysis.Metrics/)

* [Code metrics: Maintainability Index](https://learn.microsoft.com/en-us/visualstudio/code-quality/code-metrics-maintainability-index-range-and-meaning)

* [NDepend](https://www.ndepend.com/) is commercial software for performing code analysis on .NET code bases and has many advanced features that make it worth considering for organizations that wish to track code quality and who can afford the cost. The [NDepend Sample Reports](https://www.ndepend.com/sample-reports/) demonstrate useful ways to report code analysis metrics.

* This page documents findings originally discussed in [ScottPlot issue #2454](https://github.com/ScottPlot/ScottPlot/issues/2454)
November 7th, 2022

Creating Bitmaps from Scratch in C#

This project how to represent bitmap data in a plain old C object (POCO) to create images from scratch using C# and no dependencies. Common graphics libraries like SkiaSharp, ImageSharp, System.Drawing, and Maui.Graphics can read and write bitmaps in memory, so a POCO that stores image data and converts it to a bitmap byte allows creation of platform-agnostic APIs that can be interfaced from any graphics library.

This page demonstrates how to use C# (.NET 6.0) to create bitmap images from scratch. Bitmap images can then be saved to disk and viewed with any image editing program, or they can consumed as a byte array in memory by a graphics library. There are various bitmap image formats (grayscale, indexed colors, 16-bit, 32-bit, transparent, etc.) but code here demonstrates the simplest common case (8-bit RGB color).

Representing Color

The following struct represents RGB color as 3 byte values and has helper methods for creating new colors.

public struct RawColor
{
    public readonly byte R, G, B;

    public RawColor(byte r, byte g, byte b)
    {
        (R, G, B) = (r, g, b);
    }

    public static RawColor Random(Random rand)
    {
        byte r = (byte)rand.Next(256);
        byte g = (byte)rand.Next(256);
        byte b = (byte)rand.Next(256);
        return new RawColor(r, g, b);
    }

    public static RawColor Gray(byte value)
    {
        return new RawColor(value, value, value);
    }
}

A color class like this could be extended to support additional niceties. Refer to SkiaSharp's SKColor.cs, System.Drawing's Color.cs, and Maui.Graphics' Color.cs for examples and implementation details. I commonly find the following features useful include when writing a color class:

  • A static class with named colors e.g., RawColors.Blue
  • Conversion to/from ARGB e.g., RawColor.FromAGRB(123456)
  • Conversion to/from HTML e.g., RawColor.FromHtml(#003366)
  • Conversion between RGB and HSL/HSV
  • Helper functions to Lighten() and Darken()
  • Helper functions to ShiftHue()
  • Extension methods to convert to common other formats like SKColor

Representing the Bitmap Image

This is the entire image class and it serves a few specific roles:

  • Store image data in a byte array arranged identically to how it will be exported in the bitmap
  • Provide helper methods to get/set pixel color
  • Provide a method to return the image as a bitmap by adding a minimal header
public class RawBitmap
{
    public readonly int Width;
    public readonly int Height;
    private readonly byte[] ImageBytes;

    public RawBitmap(int width, int height)
    {
        Width = width;
        Height = height;
        ImageBytes = new byte[width * height * 4];
    }

    public void SetPixel(int x, int y, RawColor color)
    {
        int offset = ((Height - y - 1) * Width + x) * 4;
        ImageBytes[offset + 0] = color.B;
        ImageBytes[offset + 1] = color.G;
        ImageBytes[offset + 2] = color.R;
    }

    public byte[] GetBitmapBytes()
    {
        const int imageHeaderSize = 54;
        byte[] bmpBytes = new byte[ImageBytes.Length + imageHeaderSize];
        bmpBytes[0] = (byte)'B';
        bmpBytes[1] = (byte)'M';
        bmpBytes[14] = 40;
        Array.Copy(BitConverter.GetBytes(bmpBytes.Length), 0, bmpBytes, 2, 4);
        Array.Copy(BitConverter.GetBytes(imageHeaderSize), 0, bmpBytes, 10, 4);
        Array.Copy(BitConverter.GetBytes(Width), 0, bmpBytes, 18, 4);
        Array.Copy(BitConverter.GetBytes(Height), 0, bmpBytes, 22, 4);
        Array.Copy(BitConverter.GetBytes(32), 0, bmpBytes, 28, 2);
        Array.Copy(BitConverter.GetBytes(ImageBytes.Length), 0, bmpBytes, 34, 4);
        Array.Copy(ImageBytes, 0, bmpBytes, imageHeaderSize, ImageBytes.Length);
        return bmpBytes;
    }

    public void Save(string filename)
    {
        byte[] bytes = GetBitmapBytes();
        File.WriteAllBytes(filename, bytes);
    }
}

Generate Images from Scratch

The following code uses the bitmap class and color struct above to create test images

Random Colors

RawBitmap bmp = new(400, 300);
Random rand = new();
for (int x = 0; x < bmp.Width; x++)
    for (int y = 0; y < bmp.Height; y++)
        bmp.SetPixel(x, y, RawColor.Random(rand));
bmp.Save("random-rgb.bmp");

Rainbow

RawBitmap bmp = new(400, 300);
Random rand = new();
for (int x = 0; x < bmp.Width; x++)
{
    for (int y = 0; y < bmp.Height; y++)
    {
        byte r = (byte)(255.0 * x / bmp.Width);
        byte g = (byte)(255.0 * y / bmp.Height);
        byte b = (byte)(255 - 255.0 * x / bmp.Width);
        RawColor color = new(r, g, b);
        bmp.SetPixel(x, y, color);
    }
}
bmp.Save("rainbow.bmp");

Rectangles

RawBitmap bmp = new(400, 300);
Random rand = new();
for (int i = 0; i < 1000; i++)
{
    int rectX = rand.Next(bmp.Width);
    int rectY = rand.Next(bmp.Height);
    int rectWidth = rand.Next(50);
    int rectHeight = rand.Next(50);
    RawColor color = RawColor.Random(rand);

    for (int x = rectX; x < rectX + rectWidth; x++)
    {
        for (int y = rectY; y < rectY + rectHeight; y++)
        {
            if (x < 0 || x >= bmp.Width) continue;
            if (y < 0 || y >= bmp.Height) continue;
            bmp.SetPixel(x, y, color);
        }
    }
}
bmp.Save("rectangles.bmp");

Interfacing Graphics Libraries

The following code demonstrates how to load the bitmap byte arrays generated above into common graphics libraries and save the result as a JPEG file. Although the bitmap byte array can be written directly to disk as a .bmp file, these third-party libraries are required to encode images in additional formats like JPEG.

System.Drawing

using System.Drawing;

static void SaveBitmap(byte[] bytes, string filename = "demo.jpg")
{
    using MemoryStream ms = new(bytes);
    using Image img = Bitmap.FromStream(ms);
    img.Save(filename);
}

SkiaSharp

using SkiaSharp;

static void SaveBitmap(byte[] bytes, string filename = "demo.jpg")
{
    using SKBitmap bmp = SKBitmap.Decode(bytes);
    using SKFileWStream fs = new(filename);
    bmp.Encode(fs, SKEncodedImageFormat.Jpeg, quality: 95);
}

ImageSharp

using SixLabors.ImageSharp;
using SixLabors.ImageSharp.Formats.Jpeg;

static void SaveBitmap(byte[] bytes, string filename = "demo.jpg")
{
    using Image img = Image.Load(bytes);
    JpegEncoder encoder = new() { Quality = 95 };
    img.Save(filename, encoder);
}

Resources

Markdown source code last modified on November 8th, 2022
---
Title: Creating Bitmaps from Scratch in C#
Description: How to create a bitmap and set pixel colors in memory and save the result to disk or convert it to a traditional image format
Date: 2022-11-07 18:45PM EST
Tags: csharp, graphics
---

# Creating Bitmaps from Scratch in C# 

**This project how to represent bitmap data in a plain old C object (POCO) to create images from scratch using C# and no dependencies.** 
Common graphics libraries like [SkiaSharp](https://swharden.com/csdv/skiasharp/), [ImageSharp](https://swharden.com/csdv/platforms/imagesharp/), [System.Drawing](https://swharden.com/csdv/system.drawing/), and [Maui.Graphics](https://swharden.com/csdv/maui.graphics/) can read and write bitmaps in memory, so a [POCO](https://en.wikipedia.org/wiki/POCO) that stores image data and converts it to a bitmap byte allows creation of platform-agnostic APIs that can be interfaced from any graphics library.

**This page demonstrates how to use C# (.NET 6.0) to create bitmap images from scratch.** Bitmap images can then be saved to disk and viewed with any image editing program, or they can consumed as a byte array in memory by a graphics library. There are various [bitmap image formats](https://learn.microsoft.com/en-us/dotnet/desktop/winforms/advanced/types-of-bitmaps?view=netframeworkdesktop-4.8) (grayscale, indexed colors, 16-bit, 32-bit, transparent, etc.) but code here demonstrates the simplest common case (8-bit RGB color).

## Representing Color

The following `struct` represents RGB color as 3 `byte` values and has helper methods for creating new colors. 

```cs
public struct RawColor
{
    public readonly byte R, G, B;

    public RawColor(byte r, byte g, byte b)
    {
        (R, G, B) = (r, g, b);
    }

    public static RawColor Random(Random rand)
    {
        byte r = (byte)rand.Next(256);
        byte g = (byte)rand.Next(256);
        byte b = (byte)rand.Next(256);
        return new RawColor(r, g, b);
    }

    public static RawColor Gray(byte value)
    {
        return new RawColor(value, value, value);
    }
}
```

A color class like this could be extended to support additional niceties. Refer to [SkiaSharp's `SKColor.cs`](https://github.com/mono/SkiaSharp/blob/main/binding/Binding/SKColor.cs), [System.Drawing's `Color.cs`](https://github.com/dotnet/runtime/blob/main/src/libraries/System.Drawing.Primitives/src/System/Drawing/Color.cs), and [Maui.Graphics' `Color.cs`](https://github.com/dotnet/maui/blob/main/src/Graphics/src/Graphics/Color.cs) for examples and implementation details. I commonly find the following features useful include when writing a color class:

* A static class with named colors e.g., `RawColors.Blue`
* Conversion to/from ARGB e.g., `RawColor.FromAGRB(123456)`
* Conversion to/from HTML e.g., `RawColor.FromHtml(#003366)`
* Conversion between RGB and [HSL/HSV](https://en.wikipedia.org/wiki/HSL_and_HSV)
* Helper functions to `Lighten()` and `Darken()`
* Helper functions to `ShiftHue()`
* Extension methods to convert to common other formats like `SKColor`

## Representing the Bitmap Image

This is the entire image class and it serves a few specific roles:

* Store image data in a byte array arranged identically to how it will be exported in the bitmap
* Provide helper methods to get/set pixel color
* Provide a method to return the image as a bitmap by adding a minimal header

```cs
public class RawBitmap
{
    public readonly int Width;
    public readonly int Height;
    private readonly byte[] ImageBytes;

    public RawBitmap(int width, int height)
    {
        Width = width;
        Height = height;
        ImageBytes = new byte[width * height * 4];
    }

    public void SetPixel(int x, int y, RawColor color)
    {
        int offset = ((Height - y - 1) * Width + x) * 4;
        ImageBytes[offset + 0] = color.B;
        ImageBytes[offset + 1] = color.G;
        ImageBytes[offset + 2] = color.R;
    }

    public byte[] GetBitmapBytes()
    {
        const int imageHeaderSize = 54;
        byte[] bmpBytes = new byte[ImageBytes.Length + imageHeaderSize];
        bmpBytes[0] = (byte)'B';
        bmpBytes[1] = (byte)'M';
        bmpBytes[14] = 40;
        Array.Copy(BitConverter.GetBytes(bmpBytes.Length), 0, bmpBytes, 2, 4);
        Array.Copy(BitConverter.GetBytes(imageHeaderSize), 0, bmpBytes, 10, 4);
        Array.Copy(BitConverter.GetBytes(Width), 0, bmpBytes, 18, 4);
        Array.Copy(BitConverter.GetBytes(Height), 0, bmpBytes, 22, 4);
        Array.Copy(BitConverter.GetBytes(32), 0, bmpBytes, 28, 2);
        Array.Copy(BitConverter.GetBytes(ImageBytes.Length), 0, bmpBytes, 34, 4);
        Array.Copy(ImageBytes, 0, bmpBytes, imageHeaderSize, ImageBytes.Length);
        return bmpBytes;
    }

    public void Save(string filename)
    {
        byte[] bytes = GetBitmapBytes();
        File.WriteAllBytes(filename, bytes);
    }
}
```

## Generate Images from Scratch

The following code uses the bitmap class and color struct above to create test images

### Random Colors

```cs
RawBitmap bmp = new(400, 300);
Random rand = new();
for (int x = 0; x < bmp.Width; x++)
    for (int y = 0; y < bmp.Height; y++)
        bmp.SetPixel(x, y, RawColor.Random(rand));
bmp.Save("random-rgb.bmp");
```

![](SkiaSharp-random-rgb.jpg)

### Rainbow
```cs
RawBitmap bmp = new(400, 300);
Random rand = new();
for (int x = 0; x < bmp.Width; x++)
{
    for (int y = 0; y < bmp.Height; y++)
    {
        byte r = (byte)(255.0 * x / bmp.Width);
        byte g = (byte)(255.0 * y / bmp.Height);
        byte b = (byte)(255 - 255.0 * x / bmp.Width);
        RawColor color = new(r, g, b);
        bmp.SetPixel(x, y, color);
    }
}
bmp.Save("rainbow.bmp");
```

![](SkiaSharp-rainbow.jpg)

### Rectangles
```cs
RawBitmap bmp = new(400, 300);
Random rand = new();
for (int i = 0; i < 1000; i++)
{
    int rectX = rand.Next(bmp.Width);
    int rectY = rand.Next(bmp.Height);
    int rectWidth = rand.Next(50);
    int rectHeight = rand.Next(50);
    RawColor color = RawColor.Random(rand);

    for (int x = rectX; x < rectX + rectWidth; x++)
    {
        for (int y = rectY; y < rectY + rectHeight; y++)
        {
            if (x < 0 || x >= bmp.Width) continue;
            if (y < 0 || y >= bmp.Height) continue;
            bmp.SetPixel(x, y, color);
        }
    }
}
bmp.Save("rectangles.bmp");
```

![](SkiaSharp-rectangles.jpg)

## Interfacing Graphics Libraries

**The following code demonstrates how to load the bitmap byte arrays generated above into common graphics libraries and save the result as a JPEG file.** Although the bitmap byte array can be written directly to disk as a .bmp file, these third-party libraries are required to encode images in additional formats like JPEG.


### System.Drawing

```cs
using System.Drawing;

static void SaveBitmap(byte[] bytes, string filename = "demo.jpg")
{
    using MemoryStream ms = new(bytes);
    using Image img = Bitmap.FromStream(ms);
    img.Save(filename);
}
```

### SkiaSharp

```cs
using SkiaSharp;

static void SaveBitmap(byte[] bytes, string filename = "demo.jpg")
{
    using SKBitmap bmp = SKBitmap.Decode(bytes);
    using SKFileWStream fs = new(filename);
    bmp.Encode(fs, SKEncodedImageFormat.Jpeg, quality: 95);
}
```

### ImageSharp

```cs
using SixLabors.ImageSharp;
using SixLabors.ImageSharp.Formats.Jpeg;

static void SaveBitmap(byte[] bytes, string filename = "demo.jpg")
{
    using Image img = Image.Load(bytes);
    JpegEncoder encoder = new() { Quality = 95 };
    img.Save(filename, encoder);
}
```

# Resources

* [C# Data Visualization](https://swharden.com/csdv/) - Resources for visualizing data using C# and the .NET platform

* [SkiaSharp: `SKColor.cs`](https://github.com/mono/SkiaSharp/blob/main/binding/Binding/SKColor.cs)

* [System.Drawing: `Color.cs`](https://github.com/dotnet/runtime/blob/main/src/libraries/System.Drawing.Primitives/src/System/Drawing/Color.cs)

* [Maui.Graphics: `Color.cs`](https://github.com/dotnet/maui/blob/main/src/Graphics/src/Graphics/Color.cs)
October 16th, 2022

Experiments in PSK-31 Synthesis

PSK-31 is a narrow-bandwidth digital mode which encodes text as an audio tone that varies phase at a known rate. To learn more about this digital mode and solve a challenging programming problem, I'm going to write a PSK-31 encoder and decoder from scratch using the C# programming language. All code created for this project is open-source, available from my PSK Experiments GitHub Repository, and released under the permissive MIT license. This page documents my progress and notes things I learn along the way.

Encoding Bits as Phase Shifts

  • PSK31 messages have a continuous carrier tone

  • Symbols are represented by "symbols", each 1/31.25 seconds long

  • If a symbol changes phase from its previous symbol it is a 0, otherwise it is a 1

Amplitude Modulation Silences Phase Transitions

Although a continuous phase-shifting tone of constant amplitude can successfully transmit PSK31 data, the abrupt phase transitions will cause splatter. If you transmit this you will be heard, but those trying to communicate using adjacent frequencies will be highly disappointed.

Hard phase transitions (splatter) Soft phase transitions (cleaner)

To reduce spectral artifacts that result from abruptly changing phase, phase transitions are silenced by shaping the waveform envelope as a sine wave so it is silent at the transition. This way the maximum rate of phase shifts is a sine wave with a period of half the baud rate. This is why the opening of a PSK31 message (a series of logical 0 bits) sounds like two tones: It's the carrier sine wave with an envelope shaped like a sine wave with a period of 31.25/2 Hz. These two tones separated by approximately 15 Hz are visible in the spectrogram (waterfall).

Encoding Text as Bits

Unlike ASCII (8 bits per character) and RTTY (5 bits per character), BPSK uses Varicode (1-10 bits per character) to encode text. Consecutive zeros 00 separate characters, so character codes must not contain 00 and they must start ane end with a 1 bit. Messages are flanked by a preamble (repeated 0 bits) and a postamble (repeated 1 bits).

NUL 1010101011
SOH 1011011011
STX 1011101101
ETX 1101110111
EOT 1011101011
ENQ 1101011111
ACK 1011101111
BEL 1011111101
BS  1011111111
HT  11101111
LF  11101
VT  1101101111
FF  1011011101
CR  11111
SO  1101110101
SI  1110101011
DLE 1011110111
DC1 1011110101
DC2 1110101101
DC3 1110101111
DC4 1101011011
NAK 1101101011
SYN 1101101101
ETB 1101010111
CAN 1101111011
EM  1101111101
SUB 1110110111
ESC 1101010101
FS  1101011101
GS  1110111011
RS  1011111011
US  1101111111
SP  1
!   111111111
"   101011111
#   111110101
$   111011011
%   1011010101
&   1010111011
'   101111111
(   11111011
)   11110111
*   101101111
+   111011111
,   1110101
-   110101
.   1010111
/   110101111
0   10110111
1   10111101
2   11101101
3   11111111
4   101110111
5   101011011
6   101101011
7   110101101
8   110101011
9   110110111
:   11110101
;   110111101
<   111101101
=   1010101
>   111010111
?   1010101111
@   1010111101
A   1111101
B   11101011
C   10101101
D   10110101
E   1110111
F   11011011
G   11111101
H   101010101
I   1111111
J   111111101
K   101111101
L   11010111
M   10111011
N   11011101
O   10101011
P   11010101
Q   111011101
R   10101111
S   1101111
T   1101101
U   101010111
V   110110101
X   101011101
Y   101110101
Z   101111011
[   1010101101
\   111110111
]   111101111
^   111111011
_   1010111111
.   101101101
/   1011011111
a   1011
b   1011111
c   101111
d   101101
e   11
f   111101
g   1011011
h   101011
i   1101
j   111101011
k   10111111
l   11011
m   111011
n   1111
o   111
p   111111
q   110111111
r   10101
s   10111
t   101
u   110111
v   1111011
w   1101011
x   11011111
y   1011101
z   111010101
{   1010110111
|   110111011
}   1010110101
~   1011010111
DEL 1110110101

How to Generate a PSK Waveform

Now that we've covered the major steps of PSK31 message composition and modulation, let's go through the steps ot generate a PSK31 message in code.

Step 1: Convert a Message to Varicode

Here's the gist of how I store my varicode table in code. Note that the struct has an additional Description field which is useful for decoding and debugging.

public struct VaricodeSymbol
{
    public readonly string Symbol;
    public string BitString;
    public int[] Bits;
    public readonly string Description;

    public VaricodeSymbol(string symbol, string bitString, string? description = null)
    {
        Symbol = symbol;
        BitString = bitString;
        Bits = bitString.ToCharArray().Select(x => x == '1' ? 1 : 0).ToArray();
        Description = description ?? string.Empty;
    }
}
static VaricodeSymbol[] GetAllSymbols() => new VaricodeSymbol[]
{
    new("NUL", "1010101011", "Null character"),
    new("LF", "11101", "Line feed"),
    new("CR", "11111", "Carriage return"),
    new("SP", "1", "Space"),
    new("a", "1011"),
    new("b", "1011111"),
    new("c", "101111"),
    // etc...
};

I won't show how I do the message-to-varicode lookup, but it's trivial. Here's the final function I use to generate Varicode bits from a string:

static int[] GetVaricodeBits(string message)
{
    List<int> bits = new();

    // add a preamble of repeated zeros
    for (int i=0; i<20; i++)
        bits.Add(0);

    // encode each character of a message
    foreach (char character in message)
    {
        VaricodeSymbol symbol = Lookup(character);
        bits.AddRange(symbol.Bits);
        bits.AddRange(CharacterSeparator);
    }

    // add a postamble of repeated ones
    for (int i=0; i<20; i++)
        bits.Add(1);

    return bits.ToArray();
}

Step 2: Determine Phase Shifts

Now that we have our Varicode bits, we need to generate an array to indicate phase transitions. A transition occurs every time a bit changes value form the previous bit. This code returns phase as an array of double given the bits from a Varicode message.

public static double[] GetPhaseShifts(int[] bits, double phase1 = 0, double phase2 = Math.PI)
{
    double[] phases = new double[bits.Length];
    for (int i = 0; i < bits.Length; i++)
    {
        double previousPhase = i > 0 ? phases[i - 1] : phase1;
        double oppositePhase = previousPhase == phase1 ? phase2 : phase1;
        phases[i] = bits[i] == 1 ? previousPhase : oppositePhase;
    }
    return phases;
}

Step 3: Generate the Waveform

These constants will be used to define the shape of the waveform:

public const int SampleRate = 8000;
public const double Frequency = 1000;
public const double BaudRate = 31.25;

This minimal code generates a decipherable PSK-31 message, but it does not silence the phase transitions so it produces a lot of splatter. This function must be refined to shape the waveform such that phase transitions are silenced.

public double[] GetWaveformBPSK(double[] phases)
{
    int totalSamples = (int)(phases.Length * SampleRate / BaudRate);
    double[] wave = new double[totalSamples];
    for (int i = 0; i < wave.Length; i++)
    {
        double time = (double)i / SampleRate;
        int frame = (int)(time * BaudRate);
        double phaseShift = phases[frame];
        wave[i] = Math.Cos(2 * Math.PI * Frequency * time + phaseShift);
    }
    return wave;
}

Step 4: Generate the Waveform with Amplitude Modulation

This is the same function as above, but with extra logic for amplitude-modulating the waveform in the shape of a sine wave to silence phase transitions.

public double[] GetWaveformBPSK(double[] phases)
{
    int baudSamples = (int)(SampleRate / BaudRate);
    double samplesPerBit = SampleRate / BaudRate;
    int totalSamples = (int)(phases.Length * SampleRate / BaudRate);
    double[] wave = new double[totalSamples];

    // create the amplitude envelope sized for a single bit
    double[] envelope = new double[(int)samplesPerBit];
    for (int i = 0; i < envelope.Length; i++)
        envelope[i] = Math.Sin((i + .5) * Math.PI / envelope.Length);

    for (int i = 0; i < wave.Length; i++)
    {
        // phase modulated carrier
        double time = (double)i / SampleRate;
        int frame = (int)(time * BaudRate);
        double phaseShift = phases[frame];
        wave[i] = Math.Cos(2 * Math.PI * Frequency * time + phaseShift);

        // envelope at phase transitions
        int firstSample = (int)(frame * SampleRate / BaudRate);
        int distanceFromFrameStart = i - firstSample;
        int distanceFromFrameEnd = baudSamples - distanceFromFrameStart + 1;
        bool isFirstHalfOfFrame = distanceFromFrameStart < distanceFromFrameEnd;
        bool samePhaseAsLast = frame == 0 ? false : phases[frame - 1] == phases[frame];
        bool samePhaseAsNext = frame == phases.Length - 1 ? false : phases[frame + 1] == phases[frame];
        bool rampUp = isFirstHalfOfFrame && !samePhaseAsLast;
        bool rampDown = !isFirstHalfOfFrame && !samePhaseAsNext;

        if (rampUp)
            wave[i] *= envelope[distanceFromFrameStart];

        if (rampDown)
            wave[i] *= envelope[distanceFromFrameEnd];
    }

    return wave;
}

PSK31 Encoder Program

I wrapped the functionality above in a Windows Forms GUI that allows the user to type a message, specify frequency, baud rate, and whether or not to refine the envelope to reduce splatter, then either play or save the result. An interactive ScottPlot Chart allows the user to inspect the waveform.

Sample PSK-31 Transmissions

These audio files encode the text The Quick Brown Fox Jumped Over The Lazy Dog 1234567890 Times! in 1kHz BPSK at various baud rates.

Non-Standard Baud Rates

Let's see what PSK-3 sounds like. This mode encodes data at a rate of 3 bits per second. Note that Varicode characters may require up to ten bits, so this is pretty slow. On the other hand the side tones are closer to the carrier and the total bandwidth is much smaller. The message here has been shortened to just my callsign, AJ4VD.

Encode PSK-31 In Your Browser

After implementing the C# encoder described above I created a JavaScript version (as per Atwood's Law).

Decoding PSK-31

Considering all the steps for encoding PSK-31 transmissions are already described on this page, it doesn't require too much additional effort to create a decoder using basic software techniques. Once symbol phases are detected it's easy to work backwards: detect phase transitions (logical 0s) or repeats (logical 1s), treat consecutive zeros as a character separator, then look-up characters according to the Varicode table. The tricky bit is analyzing the source audio to generate the array of phase offsets.

The simplest way to decode PSK-31 transmissions leans on the fact that we already know the baud rate: 31.25 symbols per second, or one symbol every 256 samples at 8kHz sample rate. The audio signal can be segregated into many 256-sample bins, and processed by FFT. Once the center frequency is determined, the FFT power at this frequency can be calculated for each bin. Signal offset can be adjusted to minimize the imaginary component of the FFTs at the carrier frequency, then the real component will be strongly positive or negative, allowing phase transitions to be easily detected.

There are more advanced techniques to improve BPSK decoding, such as continuously adjusting frequency and phase alignment (synchronization). A Costas loop can help lock onto the carrier frequency while preserving its phase. Visit Simple BPSK31 Decoding with Python for an excellent demonstration of how to decode BPSK31 using these advanced techniques.

A crude C# implementation of a BPSK decoded is available on GitHub in the PSK Experiments repository

Encoding PSK-31 in Hardware

Since BPSK is just a carrier that applies periodic 180º phase-shifts, it's easy to generate in hardware by directly modulating the signal source. A good example of this is KA7OEI's PSK31 transmitter which feeds the output of an oscillator through an even or odd number of NAND gates (from a 74HC00) to produce two signals of opposite phase.

Quadrature Phase Shift Keying (QPSK)

Unlike the 0º and 180º phases of binary phase shift keying (BPSK), quadrature phase shift keying (QPSK) encodes extra data into each symbol by uses a larger number of phases. When QPSK-31 is used in amateur radio these extra bits aren't used to send messages faster but instead send them more reliably using convolutional coding and error correction. These additional features come at a cost (an extra 3 dB SNR is required), and in practice QPSK is not used as much by amateur radio operators.

QPSK encoding/decoding and convolutional encoding/decoding are outside the scope of this page, but excellent information exists on the Wikipedia: QPSK and in the US Naval Academy's EC314 Lesson 23: Digital Modulation document.

PSK-31 in 2022

After all that, it turns out PSK-31 isn't that popular anymore. These days it seems the FT-8 digital mode with WSJT-X software is orders of magnitude more popular ??

Resources

Markdown source code last modified on December 17th, 2022
---
Title: Experiments in PSK-31 Synthesis
Description: How to encode and decode PSK-31 messages using C#
Date: 2022-10-16 23:26PM EST
Tags: amateur radio, csharp
---

# Experiments in PSK-31 Synthesis

**PSK-31 is a narrow-bandwidth digital mode which encodes text as an audio tone that varies phase at a known rate.** To learn more about this digital mode and solve a challenging programming problem, I'm going to write a PSK-31 encoder and decoder from scratch using the C# programming language. All code created for this project is open-source, available from my [PSK Experiments GitHub Repository](https://github.com/swharden/psk-experiments), and released under the permissive MIT license. This page documents my progress and notes things I learn along the way.

<a href="psk-waterfall2.jpg"><img src="psk-waterfall2.jpg"></a>

## Encoding Bits as Phase Shifts

<!--
<img src="bpsk.png" class="w-75 d-block mx-auto my-3">
-->

* PSK31 messages have a continuous carrier tone

* Symbols are represented by "symbols", each 1/31.25 seconds long

* If a symbol changes phase from its previous symbol it is a `0`, otherwise it is a `1`

<a href="theory.png"><img src="theory.png" class="d-block mx-auto"></a>

## Amplitude Modulation Silences Phase Transitions

Although a continuous phase-shifting tone of constant amplitude can successfully transmit PSK31 data, the abrupt phase transitions will cause splatter. If you transmit this you will be heard, but those trying to communicate using adjacent frequencies will be highly disappointed.

<div class="text-center">

Hard phase transitions (splatter) | Soft phase transitions (cleaner)
---|---
<img src="splatter.png">|<img src="no-splatter.png">

</div>

To reduce spectral artifacts that result from abruptly changing phase, phase transitions are silenced by shaping the waveform envelope as a sine wave so it is silent at the transition. This way the maximum rate of phase shifts is a sine wave with a period of half the baud rate. This is why the opening of a PSK31 message (a series of logical `0` bits) sounds like two tones: It's the carrier sine wave with an envelope shaped like a sine wave with a period of 31.25/2 Hz. These two tones separated by approximately 15 Hz are visible in the spectrogram (waterfall).

![](modulation.png)

## Encoding Text as Bits

**Unlike ASCII (8 bits per character) and RTTY (5 bits per character), BPSK uses [Varicode](https://en.wikipedia.org/wiki/Varicode) (1-10 bits per character) to encode text.** Consecutive zeros `00` separate characters, so character codes must not contain `00` and they must start ane end with a `1` bit. Messages are flanked by a preamble (repeated `0` bits) and a postamble (repeated `1` bits).

<table cellpadding=30 style="white-space: nowrap; font-size: .8em;" class="m-0 p-0 mx-auto"><tr><td valign=top>
<TT>NUL 1010101011</TT>
<BR><TT>SOH 1011011011</TT>
<BR><TT>STX 1011101101</TT>
<BR><TT>ETX 1101110111</TT>
<BR><TT>EOT 1011101011</TT>
<BR><TT>ENQ 1101011111</TT>
<BR><TT>ACK 1011101111</TT>
<BR><TT>BEL 1011111101</TT>
<BR><TT>BS&nbsp; 1011111111</TT>
<BR><TT>HT&nbsp; 11101111</TT>
<BR><TT>LF&nbsp; 11101</TT>
<BR><TT>VT&nbsp; 1101101111</TT>
<BR><TT>FF&nbsp; 1011011101</TT>
<BR><TT>CR&nbsp; 11111</TT>
<BR><TT>SO&nbsp; 1101110101</TT>
<BR><TT>SI&nbsp; 1110101011</TT>
<BR><TT>DLE 1011110111</TT>
<BR><TT>DC1 1011110101</TT>
<BR><TT>DC2 1110101101</TT>
<BR><TT>DC3 1110101111</TT>
<BR><TT>DC4 1101011011</TT>
<BR><TT>NAK 1101101011</TT>
<BR><TT>SYN 1101101101</TT>
<BR><TT>ETB 1101010111</TT>
<BR><TT>CAN 1101111011</TT>
<BR><TT>EM&nbsp; 1101111101</TT>
<BR><TT>SUB 1110110111</TT>
<BR><TT>ESC 1101010101</TT>
<BR><TT>FS&nbsp; 1101011101</TT>
<BR><TT>GS&nbsp; 1110111011</TT>
<BR><TT>RS&nbsp; 1011111011</TT>
<BR><TT>US&nbsp; 1101111111</TT>
<BR><TT>SP&nbsp; 1</TT>
<BR><TT>!&nbsp;&nbsp; 111111111</TT>
<BR><TT>"&nbsp;&nbsp; 101011111</TT>
<BR><TT>#&nbsp;&nbsp; 111110101</TT>
<BR><TT>$&nbsp;&nbsp; 111011011</TT>
<BR><TT>%&nbsp;&nbsp; 1011010101</TT>
<BR><TT>&amp;&nbsp;&nbsp; 1010111011</TT>
<BR><TT>'&nbsp;&nbsp; 101111111</TT>
<BR><TT>(&nbsp;&nbsp; 11111011</TT>
<BR><TT>)&nbsp;&nbsp; 11110111</TT>
<BR><TT>*&nbsp;&nbsp; 101101111</TT>
</td><td valign=top width="33%" >
<TT>+&nbsp;&nbsp; 111011111</TT>
<BR><TT>,&nbsp;&nbsp; 1110101</TT>
<BR><TT>-&nbsp;&nbsp; 110101</TT>
<BR><TT>.&nbsp;&nbsp; 1010111</TT>
<BR><TT>/&nbsp;&nbsp; 110101111</TT>
<BR><TT>0&nbsp;&nbsp; 10110111</TT>
<BR><TT>1&nbsp;&nbsp; 10111101</TT>
<BR><TT>2&nbsp;&nbsp; 11101101</TT>
<BR><TT>3&nbsp;&nbsp; 11111111</TT>
<BR><TT>4&nbsp;&nbsp; 101110111</TT>
<BR><TT>5&nbsp;&nbsp; 101011011</TT>
<BR><TT>6&nbsp;&nbsp; 101101011</TT>
<BR><TT>7&nbsp;&nbsp; 110101101</TT>
<BR><TT>8&nbsp;&nbsp; 110101011</TT>
<BR><TT>9&nbsp;&nbsp; 110110111</TT>
<BR><TT>:&nbsp;&nbsp; 11110101</TT>
<BR><TT>;&nbsp;&nbsp; 110111101</TT>
<BR><TT>&lt;&nbsp;&nbsp; 111101101</TT>
<BR><TT>=&nbsp;&nbsp; 1010101</TT>
<BR><TT>>&nbsp;&nbsp; 111010111</TT>
<BR><TT>?&nbsp;&nbsp; 1010101111</TT>
<BR><TT>@&nbsp;&nbsp; 1010111101</TT>
<BR><TT>A&nbsp;&nbsp; 1111101</TT>
<BR><TT>B&nbsp;&nbsp; 11101011</TT>
<BR><TT>C&nbsp;&nbsp; 10101101</TT>
<BR><TT>D&nbsp;&nbsp; 10110101</TT>
<BR><TT>E&nbsp;&nbsp; 1110111</TT>
<BR><TT>F&nbsp;&nbsp; 11011011</TT>
<BR><TT>G&nbsp;&nbsp; 11111101</TT>
<BR><TT>H&nbsp;&nbsp; 101010101</TT>
<BR><TT>I&nbsp;&nbsp; 1111111</TT>
<BR><TT>J&nbsp;&nbsp; 111111101</TT>
<BR><TT>K&nbsp;&nbsp; 101111101</TT>
<BR><TT>L&nbsp;&nbsp; 11010111</TT>
<BR><TT>M&nbsp;&nbsp; 10111011</TT>
<BR><TT>N&nbsp;&nbsp; 11011101</TT>
<BR><TT>O&nbsp;&nbsp; 10101011</TT>
<BR><TT>P&nbsp;&nbsp; 11010101</TT>
<BR><TT>Q&nbsp;&nbsp; 111011101</TT>
<BR><TT>R&nbsp;&nbsp; 10101111</TT>
<BR><TT>S&nbsp;&nbsp; 1101111</TT>
<BR><TT>T&nbsp;&nbsp; 1101101</TT>
<BR><TT>U&nbsp;&nbsp; 101010111</TT>
</td><td valign=top width="33%" >
<TT>V&nbsp;&nbsp; 110110101</TT>
<BR><TT>X&nbsp;&nbsp; 101011101</TT>
<BR><TT>Y&nbsp;&nbsp; 101110101</TT>
<BR><TT>Z&nbsp;&nbsp; 101111011</TT>
<BR><TT>[&nbsp;&nbsp; 1010101101</TT>
<BR><TT>\&nbsp;&nbsp; 111110111</TT>
<BR><TT>]&nbsp;&nbsp; 111101111</TT>
<BR><TT>^&nbsp;&nbsp; 111111011</TT>
<BR><TT>_&nbsp;&nbsp; 1010111111</TT>
<BR><TT>.&nbsp;&nbsp; 101101101</TT>
<BR><TT>/&nbsp;&nbsp; 1011011111</TT>
<BR><TT>a&nbsp;&nbsp; 1011</TT>
<BR><TT>b&nbsp;&nbsp; 1011111</TT>
<BR><TT>c&nbsp;&nbsp; 101111</TT>
<BR><TT>d&nbsp;&nbsp; 101101</TT>
<BR><TT>e&nbsp;&nbsp; 11</TT>
<BR><TT>f&nbsp;&nbsp; 111101</TT>
<BR><TT>g&nbsp;&nbsp; 1011011</TT>
<BR><TT>h&nbsp;&nbsp; 101011</TT>
<BR><TT>i&nbsp;&nbsp; 1101</TT>
<BR><TT>j&nbsp;&nbsp; 111101011</TT>
<BR><TT>k&nbsp;&nbsp; 10111111</TT>
<BR><TT>l&nbsp;&nbsp; 11011</TT>
<BR><TT>m&nbsp;&nbsp; 111011</TT>
<BR><TT>n&nbsp;&nbsp; 1111</TT>
<BR><TT>o&nbsp;&nbsp; 111</TT>
<BR><TT>p&nbsp;&nbsp; 111111</TT>
<BR><TT>q&nbsp;&nbsp; 110111111</TT>
<BR><TT>r&nbsp;&nbsp; 10101</TT>
<BR><TT>s&nbsp;&nbsp; 10111</TT>
<BR><TT>t&nbsp;&nbsp; 101</TT>
<BR><TT>u&nbsp;&nbsp; 110111</TT>
<BR><TT>v&nbsp;&nbsp; 1111011</TT>
<BR><TT>w&nbsp;&nbsp; 1101011</TT>
<BR><TT>x&nbsp;&nbsp; 11011111</TT>
<BR><TT>y&nbsp;&nbsp; 1011101</TT>
<BR><TT>z&nbsp;&nbsp; 111010101</TT>
<BR><TT>{&nbsp;&nbsp; 1010110111</TT>
<BR><TT>|&nbsp;&nbsp; 110111011</TT>
<BR><TT>}&nbsp;&nbsp; 1010110101</TT>
<BR><TT>~&nbsp;&nbsp; 1011010111</TT>
<BR><TT>DEL 1110110101</TT>
</td></tr></table>

## How to Generate a PSK Waveform

Now that we've covered the major steps of PSK31 message composition and modulation, let's go through the steps ot generate a PSK31 message in code.

### Step 1: Convert a Message to Varicode

Here's the gist of how I store my varicode table in code. Note that the `struct` has an additional `Description` field which is useful for decoding and debugging.

```cs
public struct VaricodeSymbol
{
    public readonly string Symbol;
    public string BitString;
    public int[] Bits;
    public readonly string Description;

    public VaricodeSymbol(string symbol, string bitString, string? description = null)
    {
        Symbol = symbol;
        BitString = bitString;
        Bits = bitString.ToCharArray().Select(x => x == '1' ? 1 : 0).ToArray();
        Description = description ?? string.Empty;
    }
}
```

```cs
static VaricodeSymbol[] GetAllSymbols() => new VaricodeSymbol[]
{
    new("NUL", "1010101011", "Null character"),
    new("LF", "11101", "Line feed"),
    new("CR", "11111", "Carriage return"),
    new("SP", "1", "Space"),
    new("a", "1011"),
    new("b", "1011111"),
    new("c", "101111"),
    // etc...
};
```

I won't show how I do the message-to-varicode lookup, but it's trivial. Here's the final function I use to generate Varicode bits from a string:

```cs
static int[] GetVaricodeBits(string message)
{
    List<int> bits = new();

    // add a preamble of repeated zeros
    for (int i=0; i<20; i++)
        bits.Add(0);

    // encode each character of a message
    foreach (char character in message)
    {
        VaricodeSymbol symbol = Lookup(character);
        bits.AddRange(symbol.Bits);
        bits.AddRange(CharacterSeparator);
    }

    // add a postamble of repeated ones
    for (int i=0; i<20; i++)
        bits.Add(1);

    return bits.ToArray();
}
```

### Step 2: Determine Phase Shifts

Now that we have our Varicode bits, we need to generate an array to indicate phase transitions. A transition occurs every time a bit changes value form the previous bit. This code returns phase as an array of `double` given the bits from a Varicode message.

```cs
public static double[] GetPhaseShifts(int[] bits, double phase1 = 0, double phase2 = Math.PI)
{
    double[] phases = new double[bits.Length];
    for (int i = 0; i < bits.Length; i++)
    {
        double previousPhase = i > 0 ? phases[i - 1] : phase1;
        double oppositePhase = previousPhase == phase1 ? phase2 : phase1;
        phases[i] = bits[i] == 1 ? previousPhase : oppositePhase;
    }
    return phases;
}
```

### Step 3: Generate the Waveform

These constants will be used to define the shape of the waveform:

```cs
public const int SampleRate = 8000;
public const double Frequency = 1000;
public const double BaudRate = 31.25;
```

This minimal code generates a decipherable PSK-31 message, but it does not silence the phase transitions so it produces a lot of splatter. This function must be refined to shape the waveform such that phase transitions are silenced.

```cs
public double[] GetWaveformBPSK(double[] phases)
{
    int totalSamples = (int)(phases.Length * SampleRate / BaudRate);
    double[] wave = new double[totalSamples];
    for (int i = 0; i < wave.Length; i++)
    {
        double time = (double)i / SampleRate;
        int frame = (int)(time * BaudRate);
        double phaseShift = phases[frame];
        wave[i] = Math.Cos(2 * Math.PI * Frequency * time + phaseShift);
    }
    return wave;
}
```

### Step 4: Generate the Waveform with Amplitude Modulation

This is the same function as above, but with extra logic for amplitude-modulating the waveform in the shape of a sine wave to silence phase transitions.

```cs
public double[] GetWaveformBPSK(double[] phases)
{
    int baudSamples = (int)(SampleRate / BaudRate);
    double samplesPerBit = SampleRate / BaudRate;
    int totalSamples = (int)(phases.Length * SampleRate / BaudRate);
    double[] wave = new double[totalSamples];

    // create the amplitude envelope sized for a single bit
    double[] envelope = new double[(int)samplesPerBit];
    for (int i = 0; i < envelope.Length; i++)
        envelope[i] = Math.Sin((i + .5) * Math.PI / envelope.Length);

    for (int i = 0; i < wave.Length; i++)
    {
        // phase modulated carrier
        double time = (double)i / SampleRate;
        int frame = (int)(time * BaudRate);
        double phaseShift = phases[frame];
        wave[i] = Math.Cos(2 * Math.PI * Frequency * time + phaseShift);

        // envelope at phase transitions
        int firstSample = (int)(frame * SampleRate / BaudRate);
        int distanceFromFrameStart = i - firstSample;
        int distanceFromFrameEnd = baudSamples - distanceFromFrameStart + 1;
        bool isFirstHalfOfFrame = distanceFromFrameStart < distanceFromFrameEnd;
        bool samePhaseAsLast = frame == 0 ? false : phases[frame - 1] == phases[frame];
        bool samePhaseAsNext = frame == phases.Length - 1 ? false : phases[frame + 1] == phases[frame];
        bool rampUp = isFirstHalfOfFrame && !samePhaseAsLast;
        bool rampDown = !isFirstHalfOfFrame && !samePhaseAsNext;

        if (rampUp)
            wave[i] *= envelope[distanceFromFrameStart];

        if (rampDown)
            wave[i] *= envelope[distanceFromFrameEnd];
    }

    return wave;
}
```

## PSK31 Encoder Program

I wrapped the functionality above in a Windows Forms GUI that allows the user to type a message, specify frequency, baud rate, and whether or not to refine the envelope to reduce splatter, then either play or save the result. An interactive [ScottPlot Chart](https://scottplot.net) allows the user to inspect the waveform.

* **Download PSK31 Encoder: [PSK31-encoder.zip](PSK31-encoder.zip)**

* **PSK31 Encoder Source Code: [PSK Experiments on GitHub](https://github.com/swharden/psk-experiments)**

![](screenshot.png)

## Sample PSK-31 Transmissions

These audio files encode the text _The Quick Brown Fox Jumped Over The Lazy Dog 1234567890 Times!_ in 1kHz BPSK at various baud rates.

* PSK-31: [dog31.wav](dog31.wav)
* PSK-63: [dog63.wav](dog63.wav)
* PSK-125: [dog125.wav](dog125.wav)
* PSK-256: [dog256.wav](dog256.wav)

### Non-Standard Baud Rates

**Let's see what PSK-3 sounds like.** This mode encodes data at a rate of 3 bits per second. Note that Varicode characters may require up to ten bits, so this is pretty slow. On the other hand the side tones are closer to the carrier and the total bandwidth is much smaller. The message here has been shortened to just my callsign, AJ4VD.

* PSK-3: [psk3.wav](psk3.wav)

## Encode PSK-31 In Your Browser

After implementing the C# encoder described above I created a JavaScript version (as per <a href="https://en.wikipedia.org/wiki/Jeff_Atwood">Atwood's Law</a>).

* Try it on your phone or computer! [**Launch PskJS**](pskjs)

<a href="pskjs"><img src="pskjs.png" class="d-block mx-auto my-3"></a>

## Decoding PSK-31

Considering all the steps for _encoding_ PSK-31 transmissions are already described on this page, it doesn't require too much additional effort to create a _decoder_ using basic software techniques. Once symbol phases are detected it's easy to work backwards: detect phase transitions (logical 0s) or repeats (logical 1s), treat consecutive zeros as a character separator, then look-up characters according to the Varicode table. The tricky bit is analyzing the source audio to generate the array of phase offsets.

The simplest way to decode PSK-31 transmissions leans on the fact that we already know the baud rate: 31.25 symbols per second, or one symbol every 256 samples at 8kHz sample rate. The audio signal can be segregated into many 256-sample bins, and processed by FFT. Once the center frequency is determined, the FFT power at this frequency can be calculated for each bin. Signal offset can be adjusted to minimize the imaginary component of the FFTs at the carrier frequency, then the real component will be strongly positive or negative, allowing phase transitions to be easily detected.

<div class="row">
	<div class="col">
		<a href="psk31-receiver.html"><img src="py-fft.png" class="img-fluid"></a>
	</div>
	<div class="col">
		<a href="psk31-receiver.html"><img src="py-iq.png" class="img-fluid"></a>
	</div>
	<div class="col">
		<a href="psk31-receiver.html"><img src="py-eyediagram.png" class="img-fluid"></a>
	</div>
</div>

There are more advanced techniques to improve BPSK decoding, such as continuously adjusting frequency and phase alignment (synchronization). A [Costas loop](https://en.wikipedia.org/wiki/Costas_loop) can help lock onto the carrier frequency while preserving its phase. Visit [**Simple BPSK31 Decoding with Python**](psk31-receiver.html) for an excellent demonstration of how to decode BPSK31 using these advanced techniques.

A crude C# implementation of a BPSK decoded is available on GitHub in the [PSK Experiments](https://github.com/swharden/psk-experiments) repository

![](psk-decode.png)

## Encoding PSK-31 in Hardware

Since BPSK is just a carrier that applies periodic 180º phase-shifts, it's easy to generate in hardware by directly modulating the signal source. A good example of this is [KA7OEI's PSK31 transmitter](http://www.ka7oei.com/psk_bm_tx.html) which feeds the output of an oscillator through an even or odd number of NAND gates (from a [74HC00](https://www.mouser.com/datasheet/2/308/74HC00-105628.pdf)) to produce two signals of opposite phase.

![](pic-psk31.png)

## Quadrature Phase Shift Keying (QPSK)

**Unlike the 0º and 180º phases of binary phase shift keying (BPSK), quadrature phase shift keying (QPSK) encodes extra data into each symbol by uses a larger number of phases.** When QPSK-31 is used in amateur radio these extra bits aren't used to send messages faster but instead send them more reliably using convolutional coding and error correction. These additional features come at a cost (an extra 3 dB SNR is required), and in practice QPSK is not used as much by amateur radio operators.

QPSK encoding/decoding and convolutional encoding/decoding are outside the scope of this page, but excellent information exists on the [Wikipedia: QPSK](https://en.wikipedia.org/wiki/Phase-shift_keying#Quadrature_phase-shift_keying_(QPSK)) and in the US Naval Academy's [EC314 Lesson 23: Digital Modulation](https://www.usna.edu/ECE/ec312/Lessons/wireless/EC312_Lesson_23_Digital_Modulation_Course_Notes.pdf) document.

<a href="qpsk.png"><img src="qpsk.png" class="mx-auto d-block"></a>

## PSK-31 in 2022

After all that, it turns out PSK-31 isn't that popular anymore. These days it seems the [FT-8 digital mode](https://en.wikipedia.org/wiki/FT8) with [WSJT-X software](https://physics.princeton.edu/pulsar/k1jt/wsjtx.html) is orders of magnitude more popular ??

## Resources

* [PSK Experiments](https://github.com/swharden/psk-experiments) (GitHub) - Source code for the project shown on this page

* Software: [digipan](https://www.apkfollow.com/articles/2020/06/digipan.net.html) - A Freeware Program for PSK31 and PSK63

* Software: [fldigi](http://www.w1hkj.com/) - Supports PSK31 and other digital modes

* Software: [Digital Master 780](https://swharden.com/blog/2022-10-14-ham-radio-deluxe) - A Windows application that supports PSK31 and other digital modes bundled with Ham Radio Deluxe 5.0 (the last free version)

* Software: [WinPSK](https://www.moetronix.com/ae4jy/winpsk.htm) - open source PSK31 software for Windows

* Software: [PSKCore DLL](http://www.moetronix.com/ae4jy/pskcoredll.htm) - A Windows DLL that can be included in other software to add support for PSK31

* Software: [jacobwgillespie/psk31](https://github.com/jacobwgillespie/psk31) - Example PSK31 message generate using JavaScript

* [Digital Modulation](https://www.usna.edu/ECE/ec312/Lessons/wireless/EC312_Lesson_23_Digital_Modulation_Course_Notes.pdf) (US Naval Academy, EC314 Lesson 23) - A good description of quadrature PSK and higher order phase-shift encoding.

* [PSK-31 Specification](http://www.arrl.org/psk31-spec) (ARRL) - theory, varicode table, and convolutional code table.

* [PSK31 Description](http://aintel.bi.ehu.es/psk31.html) by G3PLX is the original / official description of the digital mode.

* [PSK31: A New Radio-Teletype Mode](http://www.arrl.org/files/file/Technology/tis/info/pdf/x9907003.pdf) (1999) by Peter Martinez, G3PLX

* [PSK31 The Easy Way](https://www.vic.wicen.org.au/wp-content/uploads/2012/05/psk31.pdf) (1999) by Alan Gibbs, VK6PG

* [Wikipedia: Varicode](https://en.wikipedia.org/wiki/Varicode) includes a table of all symbols

* [Wikipedia: QPSK](https://en.wikipedia.org/wiki/Phase-shift_keying#Quadrature_phase-shift_keying_(QPSK))

* [PSK31 Fundamentals](http://aintel.bi.ehu.es/psk31theory.html) and [PSK31 Setup](https://myplace.frontier.com/~nb6z/psk31.htm) by Peter Martinez, G3PLX

* [Varicode](http://math0.wvstateu.edu/~baker/cs240/info/varicode.html) - West Virginia State University CS240

* [Introduction to PSK31](https://sites.google.com/site/psk31matlabproject/introduction-to-psk) by engineering students at Walla Walla University

* [GNURadio PSK31 Decoder](https://sdradventure.wordpress.com/2011/10/15/gnuradio-psk31-decoder-part-1/) by VA7STH

* [Simple BPSK31](psk31-receiver.html) - a fantastic Jupyter notebook demonstrating BPSK decoding with Python

* [PySDR: Digital Modulation](https://pysdr.org/content/digital_modulation.html) - a summary of signal modulation types

* [A PIC-Based PSK31 exciter using a Balanced Modulator](http://www.ka7oei.com/psk_bm_tx.html) by Clint Turner, KA7OEI
June 23rd, 2022

Resample Time Series Data using Cubic Spline Interpolation

Cubic spline interpolation can be used to modify the sample rate of time series data. This page describes how I achieve signal resampling with spline interpolation in pure C# without any external dependencies. This technique can be used to:

  • Convert unevenly-sampled data to a series of values with a fixed sample rate

  • Convert time series data from one sample rate to another sample rate

  • Fill-in missing values from a collection of measurements

Simulating Data Samples

To simulate unevenly-sampled data I create a theoretical signal then sample it at 20 random time points.

// this function represents the signal being measured
static double f(double x) => Math.Sin(x * 10) + Math.Sin(x * 13f);

// randomly sample values from 20 time points
Random rand = new(123);
double[] sampleXs = Enumerable.Range(0, 20)
    .Select(x => rand.NextDouble())
    .OrderBy(x => x)
    .ToArray();
double[] sampleYs = sampleXs.Select(x => f(x)).ToArray();

Resample for Evenly-Spaced Data

I then generate an interpolated spline using my sampled data points as the input.

  • I can control the sample rate by defining the number of points generated in the output signal.

  • Full source code is at the bottom of this article.

(double[] xs, double[] ys) = Cubic.Interpolate1D(sampleXs, sampleYs, count: 50);

The generated points line-up perfectly with the sampled data.

There is slight deviation from the theoretical signal (and it's larger where there is more missing data) but this is an unsurprising result considering the original samples had large gaps of missing data.

Source Code

Interpolation.cs

public static class Interpolation
{
    public static (double[] xs, double[] ys) Interpolate1D(double[] xs, double[] ys, int count)
    {
        if (xs is null || ys is null || xs.Length != ys.Length)
            throw new ArgumentException($"{nameof(xs)} and {nameof(ys)} must have same length");

        int inputPointCount = xs.Length;
        double[] inputDistances = new double[inputPointCount];
        for (int i = 1; i < inputPointCount; i++)
            inputDistances[i] = inputDistances[i - 1] + xs[i] - xs[i - 1];

        double meanDistance = inputDistances.Last() / (count - 1);
        double[] evenDistances = Enumerable.Range(0, count).Select(x => x * meanDistance).ToArray();
        double[] xsOut = Interpolate(inputDistances, xs, evenDistances);
        double[] ysOut = Interpolate(inputDistances, ys, evenDistances);
        return (xsOut, ysOut);
    }

    private static double[] Interpolate(double[] xOrig, double[] yOrig, double[] xInterp)
    {
        (double[] a, double[] b) = FitMatrix(xOrig, yOrig);

        double[] yInterp = new double[xInterp.Length];
        for (int i = 0; i < yInterp.Length; i++)
        {
            int j;
            for (j = 0; j < xOrig.Length - 2; j++)
                if (xInterp[i] <= xOrig[j + 1])
                    break;

            double dx = xOrig[j + 1] - xOrig[j];
            double t = (xInterp[i] - xOrig[j]) / dx;
            double y = (1 - t) * yOrig[j] + t * yOrig[j + 1] +
                t * (1 - t) * (a[j] * (1 - t) + b[j] * t);
            yInterp[i] = y;
        }

        return yInterp;
    }

    private static (double[] a, double[] b) FitMatrix(double[] x, double[] y)
    {
        int n = x.Length;
        double[] a = new double[n - 1];
        double[] b = new double[n - 1];
        double[] r = new double[n];
        double[] A = new double[n];
        double[] B = new double[n];
        double[] C = new double[n];

        double dx1, dx2, dy1, dy2;

        dx1 = x[1] - x[0];
        C[0] = 1.0f / dx1;
        B[0] = 2.0f * C[0];
        r[0] = 3 * (y[1] - y[0]) / (dx1 * dx1);

        for (int i = 1; i < n - 1; i++)
        {
            dx1 = x[i] - x[i - 1];
            dx2 = x[i + 1] - x[i];
            A[i] = 1.0f / dx1;
            C[i] = 1.0f / dx2;
            B[i] = 2.0f * (A[i] + C[i]);
            dy1 = y[i] - y[i - 1];
            dy2 = y[i + 1] - y[i];
            r[i] = 3 * (dy1 / (dx1 * dx1) + dy2 / (dx2 * dx2));
        }

        dx1 = x[n - 1] - x[n - 2];
        dy1 = y[n - 1] - y[n - 2];
        A[n - 1] = 1.0f / dx1;
        B[n - 1] = 2.0f * A[n - 1];
        r[n - 1] = 3 * (dy1 / (dx1 * dx1));

        double[] cPrime = new double[n];
        cPrime[0] = C[0] / B[0];
        for (int i = 1; i < n; i++)
            cPrime[i] = C[i] / (B[i] - cPrime[i - 1] * A[i]);

        double[] dPrime = new double[n];
        dPrime[0] = r[0] / B[0];
        for (int i = 1; i < n; i++)
            dPrime[i] = (r[i] - dPrime[i - 1] * A[i]) / (B[i] - cPrime[i - 1] * A[i]);

        double[] k = new double[n];
        k[n - 1] = dPrime[n - 1];
        for (int i = n - 2; i >= 0; i--)
            k[i] = dPrime[i] - cPrime[i] * k[i + 1];

        for (int i = 1; i < n; i++)
        {
            dx1 = x[i] - x[i - 1];
            dy1 = y[i] - y[i - 1];
            a[i - 1] = k[i - 1] * dx1 - dy1;
            b[i - 1] = -k[i] * dx1 + dy1;
        }

        return (a, b);
    }
}

Program.cs

This is the source code I used to generate the figures on this page.

Plots were generated using ScottPlot.NET.

// this function represents the signal being measured
static double f(double x) => Math.Sin(x * 10) + Math.Sin(x * 13f);

// create points representing randomly sampled time points of a smooth curve
Random rand = new(123);
double[] sampleXs = Enumerable.Range(0, 20)
    .Select(x => rand.NextDouble())
    .OrderBy(x => x)
    .ToArray();
double[] sampleYs = sampleXs.Select(x => f(x)).ToArray();

// use 1D interpolation to create an evenly sampled curve from unevenly sampled data
(double[] xs, double[] ys) = Interpolation.Interpolate1D(sampleXs, sampleYs, count: 50);

var plt = new ScottPlot.Plot(600, 400);

double[] theoreticalXs = ScottPlot.DataGen.Range(xs.Min(), xs.Max(), .01);
double[] theoreticalYs = theoreticalXs.Select(x => f(x)).ToArray();
var perfectPlot = plt.AddScatterLines(theoreticalXs, theoreticalYs);
perfectPlot.Label = "theoretical signal";
perfectPlot.Color = plt.Palette.GetColor(2);
perfectPlot.LineStyle = ScottPlot.LineStyle.Dash;

var samplePlot = plt.AddScatterPoints(sampleXs, sampleYs);
samplePlot.Label = "sampled points";
samplePlot.Color = plt.Palette.GetColor(0);
samplePlot.MarkerSize = 10;
samplePlot.MarkerShape = ScottPlot.MarkerShape.openCircle;
samplePlot.MarkerLineWidth = 2;

var smoothPlot = plt.AddScatter(xs, ys);
smoothPlot.Label = "interpolated points";
smoothPlot.Color = plt.Palette.GetColor(3);
smoothPlot.MarkerShape = ScottPlot.MarkerShape.filledCircle;

plt.Legend();
plt.SaveFig("output.png");

2D and 3D Spline Interpolation

The interpolation method described above only considered the horizontal axis when generating evenly-spaced time points (1D interpolation). For information and code examples regarding 2D and 3D cubic spline interpolation, see my previous blog post: Spline Interpolation with C#

Resources

Markdown source code last modified on June 24th, 2022
---
Title: Resample Time Series Data using Cubic Spline Interpolation
Description: This page describes how I achieve signal resampling with spline interpolation in pure C# without any external dependencies.
Date: 2022-06-23 21:15PM EST
Tags: csharp, graphics
---

# Resample Time Series Data using Cubic Spline Interpolation

**Cubic spline interpolation can be used to modify the sample rate of time series data.** This page describes how I achieve signal resampling with spline interpolation in pure C# without any external dependencies. This technique can be used to:

* Convert unevenly-sampled data to a series of values with a fixed sample rate

* Convert time series data from one sample rate to another sample rate

* Fill-in missing values from a collection of measurements

### Simulating Data Samples

To simulate unevenly-sampled data I create a theoretical signal then sample it at 20 random time points.

```cs
// this function represents the signal being measured
static double f(double x) => Math.Sin(x * 10) + Math.Sin(x * 13f);

// randomly sample values from 20 time points
Random rand = new(123);
double[] sampleXs = Enumerable.Range(0, 20)
    .Select(x => rand.NextDouble())
    .OrderBy(x => x)
    .ToArray();
double[] sampleYs = sampleXs.Select(x => f(x)).ToArray();
```

<img src="1-samples-only.png" class="mx-auto d-block mb-5">

### Resample for Evenly-Spaced Data

I then generate an interpolated spline using my sampled data points as the input. 

* I can control the sample rate by defining the number of points generated in the output signal. 

* Full source code is at the bottom of this article.

```cs
(double[] xs, double[] ys) = Cubic.Interpolate1D(sampleXs, sampleYs, count: 50);
```

<img src="2-resample-only.png" class="mx-auto d-block mb-5">

The generated points line-up perfectly with the sampled data.

<img src="2-resample.png" class="mx-auto d-block mb-5">

There is slight deviation from the theoretical signal (and it's larger where there is more missing data) but this is an unsurprising result considering the original samples had large gaps of missing data.

<img src="3-comparison.png" class="mx-auto d-block mb-5">

## Source Code

### Interpolation.cs

```cs
public static class Interpolation
{
    public static (double[] xs, double[] ys) Interpolate1D(double[] xs, double[] ys, int count)
    {
        if (xs is null || ys is null || xs.Length != ys.Length)
            throw new ArgumentException($"{nameof(xs)} and {nameof(ys)} must have same length");

        int inputPointCount = xs.Length;
        double[] inputDistances = new double[inputPointCount];
        for (int i = 1; i < inputPointCount; i++)
            inputDistances[i] = inputDistances[i - 1] + xs[i] - xs[i - 1];

        double meanDistance = inputDistances.Last() / (count - 1);
        double[] evenDistances = Enumerable.Range(0, count).Select(x => x * meanDistance).ToArray();
        double[] xsOut = Interpolate(inputDistances, xs, evenDistances);
        double[] ysOut = Interpolate(inputDistances, ys, evenDistances);
        return (xsOut, ysOut);
    }

    private static double[] Interpolate(double[] xOrig, double[] yOrig, double[] xInterp)
    {
        (double[] a, double[] b) = FitMatrix(xOrig, yOrig);

        double[] yInterp = new double[xInterp.Length];
        for (int i = 0; i < yInterp.Length; i++)
        {
            int j;
            for (j = 0; j < xOrig.Length - 2; j++)
                if (xInterp[i] <= xOrig[j + 1])
                    break;

            double dx = xOrig[j + 1] - xOrig[j];
            double t = (xInterp[i] - xOrig[j]) / dx;
            double y = (1 - t) * yOrig[j] + t * yOrig[j + 1] +
                t * (1 - t) * (a[j] * (1 - t) + b[j] * t);
            yInterp[i] = y;
        }

        return yInterp;
    }

    private static (double[] a, double[] b) FitMatrix(double[] x, double[] y)
    {
        int n = x.Length;
        double[] a = new double[n - 1];
        double[] b = new double[n - 1];
        double[] r = new double[n];
        double[] A = new double[n];
        double[] B = new double[n];
        double[] C = new double[n];

        double dx1, dx2, dy1, dy2;

        dx1 = x[1] - x[0];
        C[0] = 1.0f / dx1;
        B[0] = 2.0f * C[0];
        r[0] = 3 * (y[1] - y[0]) / (dx1 * dx1);

        for (int i = 1; i < n - 1; i++)
        {
            dx1 = x[i] - x[i - 1];
            dx2 = x[i + 1] - x[i];
            A[i] = 1.0f / dx1;
            C[i] = 1.0f / dx2;
            B[i] = 2.0f * (A[i] + C[i]);
            dy1 = y[i] - y[i - 1];
            dy2 = y[i + 1] - y[i];
            r[i] = 3 * (dy1 / (dx1 * dx1) + dy2 / (dx2 * dx2));
        }

        dx1 = x[n - 1] - x[n - 2];
        dy1 = y[n - 1] - y[n - 2];
        A[n - 1] = 1.0f / dx1;
        B[n - 1] = 2.0f * A[n - 1];
        r[n - 1] = 3 * (dy1 / (dx1 * dx1));

        double[] cPrime = new double[n];
        cPrime[0] = C[0] / B[0];
        for (int i = 1; i < n; i++)
            cPrime[i] = C[i] / (B[i] - cPrime[i - 1] * A[i]);

        double[] dPrime = new double[n];
        dPrime[0] = r[0] / B[0];
        for (int i = 1; i < n; i++)
            dPrime[i] = (r[i] - dPrime[i - 1] * A[i]) / (B[i] - cPrime[i - 1] * A[i]);

        double[] k = new double[n];
        k[n - 1] = dPrime[n - 1];
        for (int i = n - 2; i >= 0; i--)
            k[i] = dPrime[i] - cPrime[i] * k[i + 1];

        for (int i = 1; i < n; i++)
        {
            dx1 = x[i] - x[i - 1];
            dy1 = y[i] - y[i - 1];
            a[i - 1] = k[i - 1] * dx1 - dy1;
            b[i - 1] = -k[i] * dx1 + dy1;
        }

        return (a, b);
    }
}
```

### Program.cs

This is the source code I used to generate the figures on this page. 

Plots were generated using [ScottPlot.NET](https://scottplot.net).

```cs
// this function represents the signal being measured
static double f(double x) => Math.Sin(x * 10) + Math.Sin(x * 13f);

// create points representing randomly sampled time points of a smooth curve
Random rand = new(123);
double[] sampleXs = Enumerable.Range(0, 20)
    .Select(x => rand.NextDouble())
    .OrderBy(x => x)
    .ToArray();
double[] sampleYs = sampleXs.Select(x => f(x)).ToArray();

// use 1D interpolation to create an evenly sampled curve from unevenly sampled data
(double[] xs, double[] ys) = Interpolation.Interpolate1D(sampleXs, sampleYs, count: 50);

var plt = new ScottPlot.Plot(600, 400);

double[] theoreticalXs = ScottPlot.DataGen.Range(xs.Min(), xs.Max(), .01);
double[] theoreticalYs = theoreticalXs.Select(x => f(x)).ToArray();
var perfectPlot = plt.AddScatterLines(theoreticalXs, theoreticalYs);
perfectPlot.Label = "theoretical signal";
perfectPlot.Color = plt.Palette.GetColor(2);
perfectPlot.LineStyle = ScottPlot.LineStyle.Dash;

var samplePlot = plt.AddScatterPoints(sampleXs, sampleYs);
samplePlot.Label = "sampled points";
samplePlot.Color = plt.Palette.GetColor(0);
samplePlot.MarkerSize = 10;
samplePlot.MarkerShape = ScottPlot.MarkerShape.openCircle;
samplePlot.MarkerLineWidth = 2;

var smoothPlot = plt.AddScatter(xs, ys);
smoothPlot.Label = "interpolated points";
smoothPlot.Color = plt.Palette.GetColor(3);
smoothPlot.MarkerShape = ScottPlot.MarkerShape.filledCircle;

plt.Legend();
plt.SaveFig("output.png");
```

## 2D and 3D Spline Interpolation

**The interpolation method described above only considered the horizontal axis** when generating evenly-spaced time points (1D interpolation). For information and code examples regarding 2D and 3D cubic spline interpolation, see my previous blog post: [Spline Interpolation with C#](https://swharden.com/blog/2022-01-22-spline-interpolation/) 

<a href="https://swharden.com/blog/2022-01-22-spline-interpolation/"><img src="https://swharden.com/blog/2022-01-22-spline-interpolation/screenshot.gif" class="mx-auto d-block mb-5"></a>

## Resources

* [2D and 3D Spline Interpolation with C#](https://swharden.com/blog/2022-01-22-spline-interpolation/) - Blog post from Jan, 2022

* [Spline Interpolation with ScottPlot](https://scottplot.net/cookbook/4.1/category/misc/#spline-interpolation) - Demonstrates additional types of interpolation: Bezier, Catmull-Rom, Chaikin, Cubic, etc. The project is open source under a MIT license.

* [Programmer's guide to polynomials and splines](https://wordsandbuttons.online/programmers_guide_to_polynomials_and_splines.html)
Pages