kopia lustrzana https://github.com/windytan/slowrx
complex mixdown for zoom fft
Oona Räisänen
rodzic
738917dfee
commit
2afb6d0ce8
340
src/dsp.cc
340
src/dsp.cc
|
|
@ -1,150 +1,116 @@
|
|||
#include "dsp.h"
|
||||
#include <cmath>
|
||||
|
||||
DSP::DSP() : m_fshift(0), m_sync_window(WINDOW_HANN511), m_fft_decim_ratio(4),
|
||||
m_freq_if(1900), m_window(16) {
|
||||
FM::FM(int rate, std::shared_ptr<CirBuffer<std::complex<double>>> cirbuf) : m_samplerate(rate),
|
||||
m_input_cirbuf(cirbuf), m_fshift(0), m_fft_decim_ratio(4), m_fft_len(kFFTLen),
|
||||
m_freq_if(1900) {
|
||||
|
||||
m_window[WINDOW_HANN95] = window::Hann(95);
|
||||
m_window[WINDOW_HANN127] = window::Hann(127);
|
||||
m_window[WINDOW_HANN255] = window::Hann(255);
|
||||
m_window[WINDOW_HANN511] = window::Hann(511);
|
||||
m_window[WINDOW_HANN1023] = window::Hann(1023);
|
||||
m_window[WINDOW_HANN2047] = window::Hann(2047);
|
||||
m_window[WINDOW_CHEB47] = {
|
||||
0.0004272315,0.0013212953,0.0032312239,0.0067664313,0.0127521667,0.0222058684,
|
||||
0.0363037629,0.0563165400,0.0835138389,0.1190416120,0.1637810511,0.2182020094,
|
||||
0.2822270091,0.3551233730,0.4354402894,0.5210045495,0.6089834347,0.6960162864,
|
||||
0.7784084484,0.8523735326,0.9143033652,0.9610404797,0.9901263448,1.0000000000,
|
||||
0.9901263448,0.9610404797,0.9143033652,0.8523735326,0.7784084484,0.6960162864,
|
||||
0.6089834347,0.5210045495,0.4354402894,0.3551233730,0.2822270091,0.2182020094,
|
||||
0.1637810511,0.1190416120,0.0835138389,0.0563165400,0.0363037629,0.0222058684,
|
||||
0.0127521667,0.0067664313,0.0032312239,0.0013212953,0.0004272315
|
||||
};
|
||||
m_window.push_back(window::Gauss(m_fft_len, 16, 10.0));
|
||||
m_window.push_back(window::Gauss(m_fft_len, 32, 10.0));
|
||||
m_window.push_back(window::Gauss(m_fft_len, 64, 10.0));
|
||||
m_window.push_back(window::Gauss(m_fft_len, 128, 10.0));
|
||||
m_window.push_back(window::Gauss(m_fft_len, 256, 10.0));
|
||||
|
||||
m_fft_inbuf = fftw_alloc_complex(FFT_LEN);
|
||||
m_fft_outbuf = fftw_alloc_complex(FFT_LEN);
|
||||
assert(m_fft_inbuf != NULL && m_fft_outbuf != NULL);
|
||||
// thanks arc.id.au
|
||||
m_lpf_coeffs =
|
||||
{0.000356, 0.000657, 0.000754, 0.000299, -0.001054, -0.003469, -0.006727,
|
||||
-0.010087, -0.012264, -0.011603, -0.006419, 0.004536, 0.021609, 0.043929,
|
||||
0.069344, 0.094687, 0.116329, 0.130909, 0.136054, 0.130909, 0.116329,
|
||||
0.094687, 0.069344, 0.043929, 0.021609, 0.004536, -0.006419, -0.011603,
|
||||
-0.012264, -0.010087, -0.006727, -0.003469, -0.001054, 0.000299, 0.000754,
|
||||
0.000657, 0.000356};
|
||||
|
||||
for (int i=0; i<FFT_LEN; i++) {
|
||||
m_fft_inbuf[i][0] = m_fft_inbuf[i][1] = 0.0;
|
||||
m_fft_outbuf[i][0] = m_fft_outbuf[i][1] = 0.0;
|
||||
/*m_window[WINDOW_HANN95]=window::Hann(95);
|
||||
m_window[WINDOW_HANN127]=window::Hann(127);
|
||||
m_window[WINDOW_HANN255]=window::Hann(255);
|
||||
m_window[WINDOW_HANN511]=window::Hann(511);
|
||||
m_window[WINDOW_HANN1023]=window::Hann(1023);
|
||||
m_window[WINDOW_HANN2047]=window::Hann(2047);
|
||||
m_window[WINDOW_CHEB47]={
|
||||
0.0004272315, 0.0013212953, 0.0032312239, 0.0067664313, 0.0127521667, 0.0222058684,
|
||||
0.0363037629, 0.0563165400, 0.0835138389, 0.1190416120, 0.1637810511, 0.2182020094,
|
||||
0.2822270091, 0.3551233730, 0.4354402894, 0.5210045495, 0.6089834347, 0.6960162864,
|
||||
0.7784084484, 0.8523735326, 0.9143033652, 0.9610404797, 0.9901263448, 1.0000000000,
|
||||
0.9901263448, 0.9610404797, 0.9143033652, 0.8523735326, 0.7784084484, 0.6960162864,
|
||||
0.6089834347, 0.5210045495, 0.4354402894, 0.3551233730, 0.2822270091, 0.2182020094,
|
||||
0.1637810511, 0.1190416120, 0.0835138389, 0.0563165400, 0.0363037629, 0.0222058684,
|
||||
0.0127521667, 0.0067664313, 0.0032312239, 0.0013212953, 0.0004272315
|
||||
};*/
|
||||
|
||||
m_fft_inbuf=fftw_alloc_complex(m_fft_len);
|
||||
m_fft_outbuf=fftw_alloc_complex(m_fft_len);
|
||||
assert(m_fft_inbuf!=NULL&&m_fft_outbuf!=NULL);
|
||||
|
||||
for(int i=0; i<m_fft_len; i++){
|
||||
m_fft_inbuf[i][0]=m_fft_inbuf[i][1]=0.0;
|
||||
m_fft_outbuf[i][0]=m_fft_outbuf[i][1]=0.0;
|
||||
}
|
||||
|
||||
m_mag = std::vector<double>(FFT_LEN);
|
||||
m_mag = std::vector<double>(m_fft_len);
|
||||
|
||||
m_fft_plan = fftw_plan_dft_1d(FFT_LEN, m_fft_inbuf, m_fft_outbuf, FFTW_FORWARD, FFTW_ESTIMATE);
|
||||
|
||||
m_input = std::make_shared<Input>();
|
||||
m_fft_plan = fftw_plan_dft_1d(m_fft_len, m_fft_inbuf, m_fft_outbuf, FFTW_FORWARD, FFTW_ESTIMATE);
|
||||
|
||||
}
|
||||
|
||||
DSP::~DSP() {
|
||||
FM::~FM() {
|
||||
fftw_free(m_fft_inbuf);
|
||||
fftw_free(m_fft_outbuf);
|
||||
}
|
||||
|
||||
std::shared_ptr<Input> DSP::getInput() {
|
||||
return m_input;
|
||||
}
|
||||
|
||||
// the current moment, windowed
|
||||
// will be written over buf
|
||||
// which MUST fit FFT_LEN_BIG * fftw_complex
|
||||
void DSP::calcWindowedFFT (WindowType win_type, int fft_len) {
|
||||
void FM::calcPowerSpectrum (int which_window) {
|
||||
|
||||
assert(fft_len >= (int)m_window[win_type].size() / m_fft_decim_ratio);
|
||||
for (int i = 0; i < m_fft_len; i++) {
|
||||
|
||||
for (int i=0; i<fft_len; i++)
|
||||
m_fft_inbuf[i][0] = m_fft_inbuf[i][1] = 0;
|
||||
std::complex<double> sample = m_input_cirbuf->at(i * m_fft_decim_ratio);
|
||||
|
||||
double if_phi = 0;
|
||||
for (int i = 0; i < MOMENT_LEN; i++) {
|
||||
|
||||
int win_i = i - MOMENT_LEN/2 + m_window[win_type].size()/2 ;
|
||||
|
||||
if (win_i >= 0 && win_i < (int)m_window[win_type].size()) {
|
||||
double a;
|
||||
fftw_complex mixed;
|
||||
a = m_input->m_cirbuf.at(i);
|
||||
|
||||
// mix to IF
|
||||
mixed[0] = a * cos(if_phi) - a * sin(if_phi);
|
||||
mixed[1] = a * cos(if_phi) + a * sin(if_phi);
|
||||
if_phi += 2 * M_PI * (-m_freq_if) / m_input->getSamplerate();
|
||||
|
||||
// LPF
|
||||
// TODO
|
||||
|
||||
// decimate
|
||||
if (win_i % m_fft_decim_ratio == 0) {
|
||||
m_fft_inbuf[win_i / m_fft_decim_ratio][0] = mixed[0] * m_window[win_type][win_i];
|
||||
m_fft_inbuf[win_i / m_fft_decim_ratio][1] = mixed[1] * m_window[win_type][win_i];
|
||||
}
|
||||
}
|
||||
m_fft_inbuf[i][0] = sample.real() * m_window[which_window][i];
|
||||
m_fft_inbuf[i][1] = sample.imag() * m_window[which_window][i];
|
||||
}
|
||||
|
||||
fftw_execute(m_fft_plan);
|
||||
|
||||
}
|
||||
|
||||
double DSP::calcPeakFreq (double minf, double maxf, WindowType wintype) {
|
||||
|
||||
//printf(" calcpeakfreq\n");
|
||||
int fft_len = FFT_LEN;
|
||||
|
||||
calcWindowedFFT(wintype, fft_len);
|
||||
|
||||
int peak_i = -1,peak_bin=0;
|
||||
int lobin = freq2bin(minf, fft_len);
|
||||
int hibin = freq2bin(maxf, fft_len);
|
||||
|
||||
for (int bin = lobin-1; bin <= hibin+1; bin++) {
|
||||
int i = (bin >= 0 ? bin : fft_len + bin);
|
||||
|
||||
m_mag.at(i) = complexMag(m_fft_outbuf[i]);
|
||||
if ( (bin >= lobin && bin <= hibin) &&
|
||||
(peak_i == -1 || m_mag.at(i) > m_mag.at(peak_i)) ) {
|
||||
peak_i = i;
|
||||
peak_bin = bin;
|
||||
}
|
||||
//printf("%6.2f ",m_mag.at(i));
|
||||
for (int i=0; i<m_fft_len; i++) {
|
||||
m_mag.at(i) = sqComplexMag(m_fft_outbuf[(i + m_fft_len/2) % m_fft_len]);
|
||||
}
|
||||
//printf ("%d(%d) ",peak_i,peak_bin);
|
||||
|
||||
double result = gaussianPeak(m_mag, peak_bin, true);
|
||||
|
||||
//printf(" %.3f ",result);
|
||||
|
||||
// In Hertz
|
||||
result = result / (fft_len*m_fft_decim_ratio) * m_input->getSamplerate() + m_freq_if + m_fshift;
|
||||
|
||||
|
||||
// cheb47 @ 44100 can't resolve <1700 Hz
|
||||
//if (result < 1700 && wintype == WINDOW_CHEB47)
|
||||
// result = calcPeakFreq (minf, maxf, WINDOW_HANN95);
|
||||
|
||||
//printf(" ret\n");
|
||||
//printf("%.1f Hz\n",result);
|
||||
return result;
|
||||
//printf("\n");
|
||||
|
||||
}
|
||||
|
||||
std::vector<double> DSP::calcBandPowerPerHz(const std::vector<std::vector<double>>& bands, WindowType wintype) {
|
||||
double FM::calcPeakFreq (double minf, double maxf, int which_window) {
|
||||
|
||||
int fft_len = FFT_LEN;
|
||||
calcPowerSpectrum(which_window);
|
||||
|
||||
calcWindowedFFT(wintype, fft_len);
|
||||
int peak_i = freq2bin(minf);
|
||||
|
||||
for (int i=0; i<m_fft_len; i++) {
|
||||
|
||||
if ( m_mag.at(i) > m_mag.at(peak_i)) {// && bin2freq(peak_i) >= minf && bin2freq(peak_i) <= maxf) {
|
||||
peak_i = i;
|
||||
}
|
||||
}
|
||||
|
||||
double result = gaussianPeak(m_mag, peak_i);
|
||||
|
||||
//printf("%d -> %f (%f Hz) -> %d\n",peak_i,result,bin2freq(result),freq2bin(bin2freq(result)));
|
||||
return bin2freq(result);
|
||||
|
||||
}
|
||||
|
||||
std::vector<double> FM::calcBandPowerPerHz(const std::vector<std::vector<double>>& bands, int which_window) {
|
||||
|
||||
calcPowerSpectrum(which_window);
|
||||
|
||||
Wave result;
|
||||
for (Wave band : bands) {
|
||||
double P = 0.0;
|
||||
double binwidth = 1.0 * m_input->getSamplerate() / fft_len;
|
||||
double binwidth = 1.0 * m_samplerate / m_fft_len;
|
||||
int nbins = 0;
|
||||
int lobin = freq2bin(band[0]+m_fshift, fft_len);
|
||||
int hibin = freq2bin(band[1]+m_fshift, fft_len);
|
||||
int lobin = freq2bin(band[0] + m_fshift);
|
||||
int hibin = freq2bin(band[1] + m_fshift);
|
||||
for (int bin = lobin; bin<=hibin; bin++) {
|
||||
int i = (bin >= 0 ? bin : fft_len + bin);
|
||||
P += pow(complexMag(m_fft_outbuf[i]), 2);
|
||||
int i = (bin >= 0 ? bin : m_fft_len + bin);
|
||||
P += m_mag.at(i);
|
||||
nbins++;
|
||||
}
|
||||
P = (binwidth*nbins == 0 ? 0.0 : P/(binwidth*nbins));
|
||||
|
|
@ -153,29 +119,33 @@ std::vector<double> DSP::calcBandPowerPerHz(const std::vector<std::vector<double
|
|||
return result;
|
||||
}
|
||||
|
||||
WindowType DSP::bestWindowFor(const ModeSpec& mode, double SNR) {
|
||||
WindowType WinType;
|
||||
|
||||
//double samplesInPixel = 1.0 * samplerate_ * ModeSpec[Mode].tScan / ModeSpec[Mode].ScanPixels;
|
||||
|
||||
if (SNR >= 23.0 /*&& mode.id != MODE_PD180 && mode.id != MODE_SDX*/) WinType = WINDOW_CHEB47;
|
||||
else if (SNR >= 12.0) WinType = WINDOW_HANN95;
|
||||
else if (SNR >= 8.0) WinType = WINDOW_HANN127;
|
||||
else if (SNR >= 5.0) WinType = WINDOW_HANN255;
|
||||
else if (SNR >= 4.0) WinType = WINDOW_HANN511;
|
||||
else WinType = WINDOW_HANN1023;
|
||||
|
||||
return WinType;
|
||||
std::vector<double> FM::getLastPowerSpectrum() const {
|
||||
return m_mag;
|
||||
}
|
||||
|
||||
double DSP::calcVideoSNR () {
|
||||
const double t = m_input->get_t();
|
||||
//WindowType FM::bestWindowFor(const ModeSpec& mode, double SNR) {
|
||||
// WindowType WinType;
|
||||
|
||||
//double samplesInPixel = 1.0 * samplerate_ * ModeSpec[Mode].tScan / ModeSpec[Mode].ScanPixels;
|
||||
/*
|
||||
if (SNR >= 23.0 && mode.id != MODE_PD180 && mode.id != MODE_SDX) WinType = 0;//WINDOW_CHEB47;
|
||||
else if (SNR >= 12.0) WinType = 1;//WINDOW_HANN95;
|
||||
else if (SNR >= 8.0) WinType = 2;//WINDOW_HANN127;
|
||||
else if (SNR >= 5.0) WinType = 3;//WINDOW_HANN255;
|
||||
else if (SNR >= 4.0) WinType = 4;//WINDOW_HANN511;
|
||||
else WinType = WINDOW_HANN1023;
|
||||
*/
|
||||
// return WinType;
|
||||
//}
|
||||
|
||||
double FM::calcVideoSNR () {
|
||||
/*const double t = m_input->get_t();
|
||||
|
||||
if (t >= m_next_snr_time) {
|
||||
std::vector<double> bands = calcBandPowerPerHz(
|
||||
{{FREQ_SYNC-1000,FREQ_SYNC-200},
|
||||
{FREQ_BLACK,FREQ_WHITE},
|
||||
{FREQ_WHITE+400, FREQ_WHITE+700}}
|
||||
{{1200.0-1000,1200.0-200},
|
||||
{1500.0,2300.0},
|
||||
{2300.0+400, 2300.0+700}}
|
||||
);
|
||||
double Pvideo_plus_noise = bands[1];
|
||||
double Pnoise_only = (bands[0] + bands[2]) / 2;
|
||||
|
|
@ -187,14 +157,15 @@ double DSP::calcVideoSNR () {
|
|||
m_next_snr_time = t + 50e-3;
|
||||
}
|
||||
|
||||
return m_SNR;
|
||||
return m_SNR;*/
|
||||
return 50.0;
|
||||
}
|
||||
|
||||
double DSP::calcSyncPower () {
|
||||
double FM::calcSyncPower () {
|
||||
std::vector<double> bands = calcBandPowerPerHz(
|
||||
{{FREQ_SYNC-50,FREQ_SYNC+50},
|
||||
{FREQ_BLACK,FREQ_WHITE}},
|
||||
m_sync_window
|
||||
{{1200.0-50,1200.0+50},
|
||||
{1500.0, 2300.0}},
|
||||
1
|
||||
);
|
||||
double sync;
|
||||
if (bands[1] == 0.0 || bands[0] > 4 * bands[1]) {
|
||||
|
|
@ -205,20 +176,32 @@ double DSP::calcSyncPower () {
|
|||
return sync;
|
||||
}
|
||||
|
||||
double DSP::calcVideoLevel (const ModeSpec& mode, bool is_adaptive) {
|
||||
WindowType win_type;
|
||||
double FM::calcVideoLevel (const ModeSpec& mode, int denoise_level) {
|
||||
|
||||
if (is_adaptive) win_type = bestWindowFor(mode, calcVideoSNR());
|
||||
else win_type = bestWindowFor(mode);
|
||||
std::vector<double> db_limits;
|
||||
|
||||
double freq = calcPeakFreq(FREQ_BLACK, FREQ_WHITE, win_type);
|
||||
return fclip((freq - FREQ_BLACK) / (FREQ_WHITE - FREQ_BLACK));
|
||||
if (denoise_level == 2)
|
||||
db_limits = {23.0, 12.0, 8.0, 5.0, 4.0};
|
||||
else if (denoise_level == 1)
|
||||
db_limits = {15.0, 9.0, 6.0, 4.0, 2.0};
|
||||
|
||||
int which_window = 0;
|
||||
double snr = calcVideoSNR();
|
||||
for (int i=0; i<(int)db_limits.size(); i++) {
|
||||
if (snr < db_limits.at(i)) {
|
||||
which_window = i;
|
||||
}
|
||||
}
|
||||
|
||||
double freq = calcPeakFreq(1500.0-50.0, 2300.0+50.0, which_window);
|
||||
//printf("%g\n",freq);
|
||||
return (freq - 1500.0) / (2300.0 - 1500.0);
|
||||
}
|
||||
|
||||
// return: refined peak x position (idx_peak-1 .. idx_peak+1)
|
||||
// return: refined peak/valley x position (idx_peak-1 .. idx_peak+1)
|
||||
double gaussianPeak (const Wave& signal, int idx_peak, bool wrap_around) {
|
||||
double y1,y2,y3;
|
||||
const int idx_last = signal.size()-1;
|
||||
const int sig_size = signal.size();
|
||||
|
||||
bool was_negative = false;
|
||||
if (wrap_around && idx_peak < 0) {
|
||||
|
|
@ -226,15 +209,14 @@ double gaussianPeak (const Wave& signal, int idx_peak, bool wrap_around) {
|
|||
idx_peak += signal.size();
|
||||
}
|
||||
|
||||
|
||||
if (idx_peak == 0) {
|
||||
y1 = signal.at(wrap_around ? idx_last : 1);
|
||||
y1 = signal.at(wrap_around ? sig_size-1 : 1);
|
||||
y2 = signal.at(0);
|
||||
y3 = signal.at(1);
|
||||
} else if (idx_peak == idx_last) {
|
||||
y1 = signal.at(idx_last-1);
|
||||
y2 = signal.at(idx_last);
|
||||
y3 = signal.at(wrap_around ? 0 : idx_last-1);
|
||||
} else if (idx_peak == sig_size-1) {
|
||||
y1 = signal.at(sig_size-2);
|
||||
y2 = signal.at(sig_size-1);
|
||||
y3 = signal.at(wrap_around ? 0 : sig_size-2);
|
||||
} else {
|
||||
y1 = signal.at(idx_peak - 1);
|
||||
y2 = signal.at(idx_peak);
|
||||
|
|
@ -253,10 +235,6 @@ double gaussianPeak (const Wave& signal, int idx_peak, bool wrap_around) {
|
|||
return refined;
|
||||
}
|
||||
|
||||
/*WindowType DSPworker::bestWindowFor(SSTVMode Mode) {
|
||||
return bestWindowFor(Mode, 20);
|
||||
}*/
|
||||
|
||||
namespace window {
|
||||
Wave Hann (int winlen) {
|
||||
Wave result(winlen);
|
||||
|
|
@ -282,13 +260,20 @@ namespace window {
|
|||
return result;
|
||||
}
|
||||
|
||||
Wave Gauss (int winlen) {
|
||||
Wave Gauss (int winlen, double ro2, double amp) {
|
||||
double a = amp / (sqrt(ro2) * sqrt(2 * M_PI));
|
||||
double b = 0.0;
|
||||
double c = sqrt(ro2);
|
||||
|
||||
Wave result(winlen);
|
||||
for (int i=0; i < winlen; i++)
|
||||
result[i] = 1;
|
||||
for (int i=0; i < winlen; i++) {
|
||||
double x = i-winlen/2;
|
||||
result[i] = fabs(x) < winlen/2 ? a * exp(-((x-b)*(x-b)) / (2*c*c)) : 0;
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
double sinc (double x) {
|
||||
|
|
@ -320,37 +305,42 @@ double Kernel::at(double xdist) const {
|
|||
val = (x >= 0 ? 1.0 : -1.0);
|
||||
}
|
||||
} else if (m_type == KERNEL_PULSE) {
|
||||
if (x >= -.5 && x < 1.5) {
|
||||
val = (x >= 0.0 && x < 1.0 ? 1.0 : 0.0);
|
||||
if (x >= -1 && x < 1) {
|
||||
val = (x >= -.5 && x < .5 ? 1.0 : -1.0);
|
||||
}
|
||||
}
|
||||
return val;
|
||||
}
|
||||
|
||||
double Kernel::getHalfWidth() const {
|
||||
double Kernel::getSupport() const {
|
||||
if (m_type == KERNEL_PULSE) {
|
||||
return 2.0 * m_scale;
|
||||
} else if (m_type == KERNEL_LANCZOS2) {
|
||||
return 2.0 * m_scale;
|
||||
} else if (m_type == KERNEL_LANCZOS3) {
|
||||
return 3.0 * m_scale;
|
||||
} else if (m_type == KERNEL_LANCZOS5) {
|
||||
return 5.0 * m_scale;
|
||||
}
|
||||
return 1.0 * m_scale;
|
||||
}
|
||||
|
||||
double complexMag (fftw_complex coeff) {
|
||||
return sqrt(pow(coeff[0],2) + pow(coeff[1],2));
|
||||
return sqrt(sqComplexMag(coeff));
|
||||
}
|
||||
double sqComplexMag (fftw_complex coeff) {
|
||||
return coeff[0]*coeff[0] + coeff[1]*coeff[1];
|
||||
}
|
||||
|
||||
double convolveSingle (const Wave& sig, const Kernel& kernel, double x, int border_treatment) {
|
||||
|
||||
const int signal_len = sig.size();
|
||||
const int kernel_halfwidth = std::round(kernel.getHalfWidth());
|
||||
const int kernel_support = std::ceil(kernel.getSupport());
|
||||
|
||||
double result = 0;
|
||||
|
||||
for (int i_kern=-kernel_halfwidth; i_kern<=kernel_halfwidth; i_kern++) {
|
||||
int i_src = x + i_kern;
|
||||
for (int i_kern=-kernel_support; i_kern<=kernel_support; i_kern++) {
|
||||
int i_src = std::round(x) + i_kern;
|
||||
double x_kern = i_src - x;
|
||||
double orig_value = 0;
|
||||
if (i_src >= 0 && i_src < signal_len) {
|
||||
|
|
@ -388,7 +378,7 @@ Wave upsample (const Wave& orig, double factor, int kern_type, int border_treatm
|
|||
int result_size = std::ceil(orig_size * factor);
|
||||
Wave result(result_size);
|
||||
|
||||
int halfwidth_samples = std::ceil(kern.getHalfWidth() * factor);
|
||||
int support_samples = std::ceil(kern.getSupport() * factor);
|
||||
|
||||
for (int i_src=-5; i_src<orig_size+5; i_src++) {
|
||||
double x_src = i_src + .5;
|
||||
|
|
@ -405,7 +395,7 @@ Wave upsample (const Wave& orig, double factor, int kern_type, int border_treatm
|
|||
if (val_src == 0.0)
|
||||
continue;
|
||||
|
||||
for (int i_dst=factor*i_src-halfwidth_samples-5; i_dst<=factor*i_src+halfwidth_samples+5; i_dst++) {
|
||||
for (int i_dst=factor*i_src-support_samples-5; i_dst<=factor*i_src+support_samples+5; i_dst++) {
|
||||
if (i_dst >= 0 && i_dst < (int)result_size) {
|
||||
double x_dst = (i_dst + .5) / factor;
|
||||
result.at(i_dst) += val_src * kern.at(x_dst - x_src);
|
||||
|
|
@ -529,11 +519,27 @@ std::tuple<bool,double,double> findMelody (const Wave& wave, const Melody& melod
|
|||
return { was_found, avg_fdiff, tshift };
|
||||
}
|
||||
|
||||
void DSP::set_fshift (double fshift) {
|
||||
void FM::set_fshift (double fshift) {
|
||||
m_fshift = fshift;
|
||||
}
|
||||
|
||||
int DSP::freq2bin (double freq, int fft_len) const {
|
||||
return ((freq-m_freq_if) / m_input->getSamplerate() * fft_len * m_fft_decim_ratio);
|
||||
void FM::setRate(int rate) {
|
||||
m_samplerate = rate;
|
||||
}
|
||||
|
||||
int FM::freq2bin (double freq) const {
|
||||
int bin = round((freq - m_freq_if - m_fshift) / (m_samplerate/m_fft_decim_ratio) * m_fft_len + m_fft_len/2);
|
||||
return bin;
|
||||
}
|
||||
|
||||
double FM::bin2freq (double bin) const {
|
||||
double f = (bin - m_fft_len/2) / (m_fft_len-1) * m_samplerate / m_fft_decim_ratio;
|
||||
return f + m_freq_if + m_fshift;
|
||||
}
|
||||
|
||||
std::complex<double> mixComplex(double input, double if_phi) {
|
||||
std::complex<double> output;
|
||||
output.real(input * cos(if_phi) - input * sin(if_phi));
|
||||
output.imag(input * cos(if_phi) + input * sin(if_phi));
|
||||
return output;
|
||||
}
|
||||
|
|
|
|||
53
src/dsp.h
53
src/dsp.h
|
|
@ -1,61 +1,54 @@
|
|||
#ifndef DSP_H_
|
||||
#define DSP_H_
|
||||
|
||||
#include "common.h"
|
||||
#include "input.h"
|
||||
#include <mutex>
|
||||
#include <thread>
|
||||
#include <complex>
|
||||
#include "fftw3.h"
|
||||
|
||||
// moment length only affects length of global delay, I/O interval,
|
||||
// and maximum window size.
|
||||
#define FFT_LEN 512
|
||||
|
||||
#define FREQ_MIN 500.0
|
||||
#define FREQ_MAX 3300.0
|
||||
#define FREQ_BLACK 1500.0
|
||||
#define FREQ_WHITE 2300.0
|
||||
#define FREQ_SYNC 1200.0
|
||||
#include "common.h"
|
||||
#include "input.h"
|
||||
|
||||
namespace window {
|
||||
Wave Hann (int);
|
||||
Wave Blackmann (int);
|
||||
Wave Rect (int);
|
||||
Wave Gauss (int);
|
||||
Wave Gauss (int,double,double amp=1.0);
|
||||
}
|
||||
|
||||
class Kernel {
|
||||
public:
|
||||
Kernel(int,double scale=1.0);
|
||||
double at(double) const;
|
||||
double getHalfWidth() const;
|
||||
double getSupport() const;
|
||||
private:
|
||||
int m_type;
|
||||
double m_scale;
|
||||
Wave m_precalc;
|
||||
};
|
||||
|
||||
class DSP {
|
||||
class FM {
|
||||
public:
|
||||
|
||||
DSP ();
|
||||
~DSP();
|
||||
FM (int rate, std::shared_ptr<CirBuffer<std::complex<double>>> cirbuf);
|
||||
~FM();
|
||||
|
||||
double calcPeakFreq (double minf, double maxf, WindowType win_type);
|
||||
std::vector<double> calcBandPowerPerHz (const std::vector<std::vector<double>>&, WindowType wintype=WINDOW_HANN2047);
|
||||
double calcPeakFreq (double minf, double maxf, int which_window=0);
|
||||
std::vector<double> calcBandPowerPerHz (const std::vector<std::vector<double>>&, int which_window=0);
|
||||
double calcVideoSNR ();
|
||||
double calcVideoLevel (const ModeSpec&, bool is_adaptive=false);
|
||||
double calcVideoLevel (const ModeSpec&, int denoise_level=0);
|
||||
double calcSyncPower ();
|
||||
void calcWindowedFFT (WindowType win_type, int fft_len);
|
||||
void calcPowerSpectrum (int which_window=0);
|
||||
std::vector<double> getLastPowerSpectrum() const;
|
||||
void setRate(int rate);
|
||||
|
||||
int freq2bin (double freq, int fft_len) const;
|
||||
WindowType bestWindowFor (const ModeSpec&, double SNR=99);
|
||||
int freq2bin (double freq) const;
|
||||
double bin2freq(double bin) const;
|
||||
//WindowType bestWindowFor (const ModeSpec&, double SNR=99);
|
||||
void set_fshift (double);
|
||||
|
||||
std::shared_ptr<Input> getInput ();
|
||||
|
||||
private:
|
||||
|
||||
int m_samplerate;
|
||||
fftw_complex* m_fft_inbuf;
|
||||
fftw_complex* m_fft_outbuf;
|
||||
fftw_plan m_fft_plan;
|
||||
|
|
@ -63,13 +56,15 @@ class DSP {
|
|||
double m_fshift;
|
||||
double m_next_snr_time;
|
||||
double m_SNR;
|
||||
WindowType m_sync_window;
|
||||
//WindowType m_sync_window;
|
||||
int m_fft_decim_ratio;
|
||||
double m_freq_if;
|
||||
const double m_freq_if;
|
||||
const int m_fft_len;
|
||||
std::vector<double> m_lpf_coeffs;
|
||||
|
||||
std::vector<Wave> m_window;
|
||||
std::shared_ptr<CirBuffer<std::complex<double>>> m_input_cirbuf;
|
||||
|
||||
std::shared_ptr<Input> m_input;
|
||||
};
|
||||
|
||||
Wave convolve (const Wave&, const Kernel&, int border_treatment=BORDER_ZERO);
|
||||
|
|
@ -82,5 +77,7 @@ Wave upsample (const Wave& orig, double factor, int kern_type, int border_
|
|||
double gaussianPeak (const Wave& signal, int idx_peak, bool wrap_around=false);
|
||||
double power (fftw_complex coeff);
|
||||
double complexMag (fftw_complex coeff);
|
||||
double sqComplexMag (fftw_complex coeff);
|
||||
std::complex<double> mixComplex(double input, double if_phi);
|
||||
|
||||
#endif
|
||||
|
|
|
|||
Ładowanie…
Reference in New Issue