144 lines
3.2 KiB
C++
144 lines
3.2 KiB
C++
/*
|
||
* © Copyright 2014 – Urheberrechtshinweis
|
||
* Alle Rechte vorbehalten / All Rights Reserved
|
||
*
|
||
* Programmcode ist urheberrechtlich geschuetzt.
|
||
* Das Urheberrecht liegt, soweit nicht ausdruecklich anders gekennzeichnet, bei Frank Ebner.
|
||
* Keine Verwendung ohne explizite Genehmigung.
|
||
* (vgl. § 106 ff UrhG / § 97 UrhG)
|
||
*/
|
||
|
||
#ifndef FIRREALFACTORY_H
|
||
#define FIRREALFACTORY_H
|
||
|
||
#include "Real.h"
|
||
|
||
namespace FIR {
|
||
|
||
namespace Real {
|
||
|
||
class Factory {
|
||
|
||
Kernel kernel;
|
||
|
||
float sRate_hz;
|
||
|
||
public:
|
||
|
||
/** ctor */
|
||
Factory(float sRate_hz) : sRate_hz(sRate_hz) {
|
||
;
|
||
}
|
||
|
||
/** frequency shift the kernel by multiplying with a frequency */
|
||
void shift(const float shift_hz) {
|
||
for (size_t i = 0; i < kernel.size(); ++i) {
|
||
const float t = (float) i / (float) sRate_hz;
|
||
const float real = std::sin(t * 2 * M_PI * shift_hz);
|
||
kernel[i] = kernel[i] * real;
|
||
}
|
||
}
|
||
|
||
/** create a lowpass filte kernel */
|
||
void lowpass(const float cutOff_hz, const int n = 50) {
|
||
|
||
kernel.clear();
|
||
|
||
for (int i = -n; i <= +n; ++i) {
|
||
const double t = (double) i / (double) sRate_hz;
|
||
const double tmp = 2 * M_PI * cutOff_hz * t;
|
||
const double val = (tmp == 0) ? (1) : (std::sin(tmp) / tmp);
|
||
const double res = val;// * 0.5f; // why 0.5?
|
||
if (res != res) {throw std::runtime_error("detected NaN");}
|
||
kernel.push_back(res);
|
||
}
|
||
|
||
}
|
||
|
||
/** apply hamming window to the filter */
|
||
void applyWindowHamming() {
|
||
const int n = (kernel.size()-1)/2;
|
||
int i = -n;
|
||
for (float& f : kernel) {
|
||
f *= winHamming(i+n, n*2);
|
||
++i;
|
||
}
|
||
}
|
||
|
||
// https://dsp.stackexchange.com/questions/4693/fir-filter-gain
|
||
/** normalize using the DC-part of the kernel */
|
||
void normalizeDC() {
|
||
float sum = 0;
|
||
for (auto f : kernel) {sum += f;}
|
||
for (auto& f : kernel) {f /= sum;}
|
||
}
|
||
|
||
// https://dsp.stackexchange.com/questions/4693/fir-filter-gain
|
||
void normalizeAC(const float freq_hz) {
|
||
|
||
const int n = (kernel.size()-1)/2;
|
||
int i = -n;
|
||
float sum = 0;
|
||
|
||
for (float f : kernel) {
|
||
const double t = (double) i / (double) sRate_hz;
|
||
const double r = 2 * M_PI * freq_hz * t;
|
||
const double s = std::sin(r);
|
||
sum += f * s;
|
||
++i;
|
||
}
|
||
|
||
for (auto& f : kernel) {f /= sum;}
|
||
|
||
}
|
||
|
||
|
||
|
||
|
||
|
||
Kernel getLowpass(const float cutOff_hz, const int n) {
|
||
lowpass(cutOff_hz, n);
|
||
applyWindowHamming();
|
||
normalizeDC();
|
||
return kernel;
|
||
}
|
||
|
||
Kernel getBandpass(const float width_hz, const float center_hz, const int n) {
|
||
lowpass(width_hz/2, n);
|
||
applyWindowHamming();
|
||
//normalizeDC();
|
||
shift(center_hz);
|
||
normalizeAC(center_hz);
|
||
return kernel;
|
||
}
|
||
|
||
void dumpKernel(const std::string& file, const std::string& varName) {
|
||
|
||
std::ofstream out(file);
|
||
out << "# name: " << varName << "\n";
|
||
out << "# type: matrix\n";
|
||
out << "# rows: " << kernel.size() << "\n";
|
||
out << "# columns: 1\n";
|
||
|
||
for (const float f : kernel) {
|
||
out << f << "\n";
|
||
}
|
||
out.close();
|
||
|
||
}
|
||
|
||
private:
|
||
|
||
/** get a value from the hamming window */
|
||
static double winHamming(const double t, const double size) {
|
||
return 0.54 - 0.46 * std::cos(2 * M_PI * t / size);
|
||
}
|
||
|
||
};
|
||
|
||
}
|
||
|
||
}
|
||
|
||
#endif // FIRREALFACTORY_H
|