disorders/multipath.cc at master · aicodix/disorders · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
/*
Multipath propagation

Copyright 2020 Ahmet Inan <inan@aicodix.de>
*/

#include <iostream>
#include <fstream>
#include <sstream>
#include <cmath>
#include <vector>
#include "complex.hh"
#include "decibel.hh"
#include "hilbert.hh"
#include "phasor.hh"
#include "wav.hh"

int main(int argc, char **argv)
{
	if (argc != 5) {
		std::cerr << "usage: " << argv[0] << " OUTPUT INPUT TABLE FACTOR" << std::endl;
		return 1;
	}

	typedef float value;
	typedef DSP::Complex<value> cmplx;

	const char *out_name = argv[1];
	if (out_name[0] == '-' && out_name[1] == 0)
		out_name = "/dev/stdout";
	const char *inp_name = argv[2];
	if (inp_name[0] == '-' && inp_name[1] == 0)
		inp_name = "/dev/stdin";
	const char *table_name = argv[3];
	value factor = DSP::idecibel(std::atof(argv[4]));

	DSP::ReadWAV<value> inp_file(inp_name);

	if (inp_file.channels() < 1 || inp_file.channels() > 2) {
		std::cerr << "Only real or analytic signal (one or two channels) supported." << std::endl;
		return 1;
	}

	DSP::WriteWAV<value> out_file(out_name, inp_file.rate(), inp_file.bits(), inp_file.channels());

	bool real = inp_file.channels() == 1;
	DSP::Hilbert<cmplx, 513> hilbert;

	typedef struct {
		cmplx camp;
		int delay;
	} Path;
	std::vector<Path> paths;

	std::ifstream table_file(table_name);

	std::string line;
	int dmax = 0;
	value power = 0;
	while (getline(table_file, line)) {
		std::istringstream iss(line);
		value ampl, msec, rad;
		iss >> ampl >> msec >> rad;
		cmplx camp = DSP::polar(ampl, rad);
		power += norm(camp);
		int delay = nearbyint(inp_file.rate() * msec / 1000);
		dmax = std::max(dmax, delay);
		paths.emplace_back(Path{camp, delay});
	}
	value ampl0 = sqrt(factor * power);
	power += ampl0 * ampl0;

	ampl0 /= sqrt(power);
	for (auto &path: paths)
		path.camp /= sqrt(power);

	cmplx *buf = new cmplx[dmax+1];

	while (out_file.good() && inp_file.good()) {
		cmplx input;
		inp_file.read(reinterpret_cast<value *>(&input), 1);
		if (real)
			input = hilbert(input.real());
		for (int i = dmax; i; --i)
			buf[i] = buf[i-1];
		buf[0] = input;
		cmplx sum = ampl0 * buf[0];
		for (const auto &path: paths)
			sum += path.camp * buf[path.delay];
		out_file.write(reinterpret_cast<value *>(&sum), 1);
	}

	return 0;
}