Test bench for doing the Short-time Fourier transform using channelizer techniques

gnuplot.txt

plot "<./testbench 3 1 1 0 | tee data.txt" u 1:2 w l t "-1", "data.txt" u 1:3 w l t "0", "data.txt" u 1:4 w l t "+1"

Makefile

CXXFLAGS = -std=c++11 -W -Wall -Ofast -fno-exceptions -fno-rtti -march=native -I../dsp
CXX = clang++ -stdlib=libc++
#CXX = g++

testbench: testbench.cc
	$(CXX) $(CXXFLAGS) $< -o $@

.PHONY: clean

clean:
	rm -f testbench

testbench.cc

/*
Test bench for doing the Short-time Fourier transform using channelizer techniques

Copyright 2021 Ahmet Inan <xdsopl@gmail.com>
*/

#include <cmath>
#include <iostream>
#include "decibel.hh"
#include "complex.hh"
#include "phasor.hh"
#include "window.hh"
#include "coeffs.hh"
#include "filter.hh"
#include "const.hh"
#include "fft.hh"

int main(int argc, char **argv)
{
	(void)argc;
	const int BINS = 1000;
	typedef float value;
	typedef DSP::Complex<value> cmplx;
	DSP::FastFourierTransform<BINS, cmplx, -1> fwd;
	//DSP::FastFourierTransform<BINS, cmplx, 1> bwd;

	DSP::Hann<value> window;
	//DSP::Hamming<value> window;
	//DSP::Lanczos<value> window;
	//DSP::Blackman<value> window;
	//DSP::Gauss<value> window(0.5);
	//DSP::Kaiser<value> window(2);
	//DSP::Rect<value> window;

	//DSP::Lanczos<value> filter;
	DSP::LowPass2<value> filter(1, BINS);

	const int FOLD = atoi(argv[1]);
	value win[FOLD*BINS];
	for (int i = 0; i < FOLD*BINS; ++i)
		win[i] = 1;

	if (atoi(argv[2]))
		for (int i = 0; i < FOLD*BINS; ++i)
			win[i] *= window(i, FOLD*BINS);
	if (atoi(argv[3]))
		for (int i = 0; i < FOLD*BINS; ++i)
			win[i] *= filter(i, FOLD*BINS);
	if (1) {
		value sum = 0;
		for (int i = 0; i < FOLD*BINS; ++i)
			sum += win[i];
		value mag = abs(sum);
		for (int i = 0; i < FOLD*BINS; ++i)
			win[i] /= mag;
	}
	if (atoi(argv[4])) {
		for (int i = 0; i < FOLD*BINS; ++i)
			std::cout << i << " " << win[i] << std::endl;
		return 0;
	}

	const int OVER = 20;
	for (int freq = -10*OVER; freq <= 10*OVER; ++freq) {
		cmplx inp[BINS], out[BINS];
		DSP::Phasor<cmplx> osc;
		osc.omega(freq, OVER*BINS);
		for (int i = 0; i < FOLD*BINS; ++i)
			inp[i%BINS] += win[i] * osc();
		fwd(out, inp);
		std::cout << freq / value(OVER);
		for (int i = -1; i <= 1; ++i)
			std::cout << " " << DSP::decibel(abs(out[(i+BINS)%BINS]));
		std::cout << std::endl;
	}
	return 0;
}