mirror of
				https://github.com/f4exb/sdrangel.git
				synced 2025-11-04 05:30:32 -05:00 
			
		
		
		
	
		
			
				
	
	
		
			343 lines
		
	
	
		
			10 KiB
		
	
	
	
		
			C++
		
	
	
	
	
	
			
		
		
	
	
			343 lines
		
	
	
		
			10 KiB
		
	
	
	
		
			C++
		
	
	
	
	
	
///////////////////////////////////////////////////////////////////////////////////
 | 
						|
// Copyright (C) 2023 Edouard Griffiths, F4EXB.                                  //
 | 
						|
//                                                                               //
 | 
						|
// This is the code from ft8mon: https://github.com/rtmrtmrtmrtm/ft8mon          //
 | 
						|
// written by Robert Morris, AB1HL                                               //
 | 
						|
// reformatted and adapted to Qt and SDRangel context                            //
 | 
						|
//                                                                               //
 | 
						|
// This program is free software; you can redistribute it and/or modify          //
 | 
						|
// it under the terms of the GNU General Public License as published by          //
 | 
						|
// the Free Software Foundation as version 3 of the License, or                  //
 | 
						|
// (at your option) any later version.                                           //
 | 
						|
//                                                                               //
 | 
						|
// This program is distributed in the hope that it will be useful,               //
 | 
						|
// but WITHOUT ANY WARRANTY; without even the implied warranty of                //
 | 
						|
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the                  //
 | 
						|
// GNU General Public License V3 for more details.                               //
 | 
						|
//                                                                               //
 | 
						|
// You should have received a copy of the GNU General Public License             //
 | 
						|
// along with this program. If not, see <http://www.gnu.org/licenses/>.          //
 | 
						|
///////////////////////////////////////////////////////////////////////////////////
 | 
						|
#include <math.h>
 | 
						|
#include <complex>
 | 
						|
#include <string>
 | 
						|
#include <algorithm>
 | 
						|
 | 
						|
#include "util/timeutil.h"
 | 
						|
#include "util.h"
 | 
						|
 | 
						|
namespace FT8 {
 | 
						|
 | 
						|
double now()
 | 
						|
{
 | 
						|
    return TimeUtil::nowus() / 1000000.0;
 | 
						|
}
 | 
						|
 | 
						|
//
 | 
						|
// Goertzel Algorithm for a Non-integer Frequency Index, Rick Lyons
 | 
						|
// https://www.dsprelated.com/showarticle/495.php
 | 
						|
//
 | 
						|
std::complex<float> goertzel(std::vector<float> v, int rate, int i0, int n, float hz)
 | 
						|
{
 | 
						|
    // float radians_per_sample = (hz * 2 * M_PI) / rate;
 | 
						|
    // float k = radians_per_sample * n;
 | 
						|
    float bin_hz = rate / (float)n;
 | 
						|
    float k = hz / bin_hz;
 | 
						|
 | 
						|
    float alpha = 2 * M_PI * k / n;
 | 
						|
    float beta = 2 * M_PI * k * (n - 1.0) / n;
 | 
						|
 | 
						|
    float two_cos_alpha = 2 * cos(alpha);
 | 
						|
    float a = cos(beta);
 | 
						|
    float b = -sin(beta);
 | 
						|
    float c = sin(alpha) * sin(beta) - cos(alpha) * cos(beta);
 | 
						|
    float d = sin(2 * M_PI * k);
 | 
						|
 | 
						|
    float w1 = 0;
 | 
						|
    float w2 = 0;
 | 
						|
 | 
						|
    for (int i = 0; i < n; i++)
 | 
						|
    {
 | 
						|
        float w0 = v[i0 + i] + two_cos_alpha * w1 - w2;
 | 
						|
        w2 = w1;
 | 
						|
        w1 = w0;
 | 
						|
    }
 | 
						|
 | 
						|
    float re = w1 * a + w2 * c;
 | 
						|
    float im = w1 * b + w2 * d;
 | 
						|
 | 
						|
    return std::complex<float>(re, im);
 | 
						|
}
 | 
						|
 | 
						|
float vmax(const std::vector<float> &v)
 | 
						|
{
 | 
						|
    float mx = 0;
 | 
						|
    int got = 0;
 | 
						|
    for (int i = 0; i < (int)v.size(); i++)
 | 
						|
    {
 | 
						|
        if (got == 0 || v[i] > mx)
 | 
						|
        {
 | 
						|
            got = 1;
 | 
						|
            mx = v[i];
 | 
						|
        }
 | 
						|
    }
 | 
						|
    return mx;
 | 
						|
}
 | 
						|
 | 
						|
std::vector<float> vreal(const std::vector<std::complex<float>> &a)
 | 
						|
{
 | 
						|
    std::vector<float> b(a.size());
 | 
						|
    for (int i = 0; i < (int)a.size(); i++)
 | 
						|
    {
 | 
						|
        b[i] = a[i].real();
 | 
						|
    }
 | 
						|
    return b;
 | 
						|
}
 | 
						|
 | 
						|
std::vector<float> vimag(const std::vector<std::complex<float>> &a)
 | 
						|
{
 | 
						|
    std::vector<float> b(a.size());
 | 
						|
    for (int i = 0; i < (int)a.size(); i++)
 | 
						|
    {
 | 
						|
        b[i] = a[i].imag();
 | 
						|
    }
 | 
						|
    return b;
 | 
						|
}
 | 
						|
 | 
						|
// generate 8-FSK, at 25 hz, bin size 6.25 hz,
 | 
						|
// 200 samples/second, 32 samples/symbol.
 | 
						|
// used as reference to detect pairs of symbols.
 | 
						|
// superseded by gfsk().
 | 
						|
std::vector<std::complex<float>> fsk_c(const std::vector<int> &syms)
 | 
						|
{
 | 
						|
    int n = syms.size();
 | 
						|
    std::vector<std::complex<float>> v(n * 32);
 | 
						|
    float theta = 0;
 | 
						|
    for (int si = 0; si < n; si++)
 | 
						|
    {
 | 
						|
        float hz = 25 + syms[si] * 6.25;
 | 
						|
        for (int i = 0; i < 32; i++)
 | 
						|
        {
 | 
						|
            v[si * 32 + i] = std::complex<float>(cos(theta), sin(theta));
 | 
						|
            theta += 2 * M_PI / (200 / hz);
 | 
						|
        }
 | 
						|
    }
 | 
						|
    return v;
 | 
						|
}
 | 
						|
 | 
						|
// copied from wsjt-x ft2/gfsk_pulse.f90.
 | 
						|
// b is 1.0 for FT4; 2.0 for FT8.
 | 
						|
float gfsk_point(float b, float t)
 | 
						|
{
 | 
						|
    float c = M_PI * sqrt(2.0 / log(2.0));
 | 
						|
    float x = 0.5 * (erf(c * b * (t + 0.5)) - erf(c * b * (t - 0.5)));
 | 
						|
    return x;
 | 
						|
}
 | 
						|
 | 
						|
// the smoothing window for gfsk.
 | 
						|
// run the window over impulses of symbol frequencies,
 | 
						|
// each impulse at the center of its symbol time.
 | 
						|
// three symbols wide.
 | 
						|
// most of the pulse is in the center symbol.
 | 
						|
// b is 1.0 for FT4; 2.0 for FT8.
 | 
						|
std::vector<float> gfsk_window(int samples_per_symbol, float b)
 | 
						|
{
 | 
						|
    std::vector<float> v(3 * samples_per_symbol);
 | 
						|
    float sum = 0;
 | 
						|
    for (int i = 0; i < (int)v.size(); i++)
 | 
						|
    {
 | 
						|
        float x = i / (float)samples_per_symbol;
 | 
						|
        x -= 1.5;
 | 
						|
        float y = gfsk_point(b, x);
 | 
						|
        v[i] = y;
 | 
						|
        sum += y;
 | 
						|
    }
 | 
						|
 | 
						|
    for (int i = 0; i < (int)v.size(); i++)
 | 
						|
    {
 | 
						|
        v[i] /= sum;
 | 
						|
    }
 | 
						|
 | 
						|
    return v;
 | 
						|
}
 | 
						|
 | 
						|
// gaussian-smoothed fsk.
 | 
						|
// the gaussian smooths the instantaneous frequencies,
 | 
						|
// so that the transitions between symbols don't
 | 
						|
// cause clicks.
 | 
						|
// gwin is gfsk_window(32, 2.0)
 | 
						|
std::vector<std::complex<float>> gfsk_c(
 | 
						|
    const std::vector<int> &symbols,
 | 
						|
    float hz0, float hz1,
 | 
						|
    float spacing, int rate, int symsamples,
 | 
						|
    float phase0,
 | 
						|
    const std::vector<float> &gwin
 | 
						|
)
 | 
						|
{
 | 
						|
    if (!((gwin.size() % 2) == 0))
 | 
						|
    {
 | 
						|
        std::vector<std::complex<float>> v(symsamples * symbols.size());
 | 
						|
        return v;
 | 
						|
    }
 | 
						|
 | 
						|
    // compute frequency for each symbol.
 | 
						|
    // generate a spike in the middle of each symbol time;
 | 
						|
    // the gaussian filter will turn it into a waveform.
 | 
						|
    std::vector<float> hzv(symsamples * (symbols.size() + 2), 0.0);
 | 
						|
    for (int bi = 0; bi < (int)symbols.size(); bi++)
 | 
						|
    {
 | 
						|
        float base_hz = hz0 + (hz1 - hz0) * (bi / (float)symbols.size());
 | 
						|
        float fr = base_hz + (symbols[bi] * spacing);
 | 
						|
        int mid = symsamples * (bi + 1) + symsamples / 2;
 | 
						|
        // the window has even size, so split the impulse over
 | 
						|
        // the two middle samples to be symmetric.
 | 
						|
        hzv[mid] = fr * symsamples / 2.0;
 | 
						|
        hzv[mid - 1] = fr * symsamples / 2.0;
 | 
						|
    }
 | 
						|
 | 
						|
    // repeat first and last symbols
 | 
						|
    for (int i = 0; i < symsamples; i++)
 | 
						|
    {
 | 
						|
        hzv[i] = hzv[i + symsamples];
 | 
						|
        hzv[symsamples * (symbols.size() + 1) + i] = hzv[symsamples * symbols.size() + i];
 | 
						|
    }
 | 
						|
 | 
						|
    // run the per-sample frequency vector through
 | 
						|
    // the gaussian filter.
 | 
						|
    int half = gwin.size() / 2;
 | 
						|
    std::vector<float> o(hzv.size());
 | 
						|
    for (int i = 0; i < (int)o.size(); i++)
 | 
						|
    {
 | 
						|
        float sum = 0;
 | 
						|
        for (int j = 0; j < (int)gwin.size(); j++)
 | 
						|
        {
 | 
						|
            int k = i - half + j;
 | 
						|
            if (k >= 0 && k < (int)hzv.size())
 | 
						|
            {
 | 
						|
                sum += hzv[k] * gwin[j];
 | 
						|
            }
 | 
						|
        }
 | 
						|
        o[i] = sum;
 | 
						|
    }
 | 
						|
 | 
						|
    // drop repeated first and last symbols
 | 
						|
    std::vector<float> oo(symsamples * symbols.size());
 | 
						|
    for (int i = 0; i < (int)oo.size(); i++)
 | 
						|
    {
 | 
						|
        oo[i] = o[i + symsamples];
 | 
						|
    }
 | 
						|
 | 
						|
    // now oo[i] contains the frequency for the i'th sample.
 | 
						|
 | 
						|
    std::vector<std::complex<float>> v(symsamples * symbols.size());
 | 
						|
    float theta = phase0;
 | 
						|
    for (int i = 0; i < (int)v.size(); i++)
 | 
						|
    {
 | 
						|
        v[i] = std::complex<float>(cos(theta), sin(theta));
 | 
						|
        float hz = oo[i];
 | 
						|
        theta += 2 * M_PI / (rate / hz);
 | 
						|
    }
 | 
						|
 | 
						|
    return v;
 | 
						|
}
 | 
						|
 | 
						|
// gaussian-smoothed fsk.
 | 
						|
// the gaussian smooths the instantaneous frequencies,
 | 
						|
// so that the transitions between symbols don't
 | 
						|
// cause clicks.
 | 
						|
// gwin is gfsk_window(32, 2.0)
 | 
						|
std::vector<float> gfsk_r(
 | 
						|
    const std::vector<int> &symbols,
 | 
						|
    float hz0, float hz1,
 | 
						|
    float spacing, int rate, int symsamples,
 | 
						|
    float phase0,
 | 
						|
    const std::vector<float> &gwin
 | 
						|
)
 | 
						|
{
 | 
						|
    if (!((gwin.size() % 2) == 0))
 | 
						|
    {
 | 
						|
        std::vector<float> v(symsamples * symbols.size());
 | 
						|
        return v;
 | 
						|
    }
 | 
						|
 | 
						|
    // compute frequency for each symbol.
 | 
						|
    // generate a spike in the middle of each symbol time;
 | 
						|
    // the gaussian filter will turn it into a waveform.
 | 
						|
    std::vector<float> hzv(symsamples * (symbols.size() + 2), 0.0);
 | 
						|
    for (int bi = 0; bi < (int)symbols.size(); bi++)
 | 
						|
    {
 | 
						|
        float base_hz = hz0 + (hz1 - hz0) * (bi / (float)symbols.size());
 | 
						|
        float fr = base_hz + (symbols[bi] * spacing);
 | 
						|
        int mid = symsamples * (bi + 1) + symsamples / 2;
 | 
						|
        // the window has even size, so split the impulse over
 | 
						|
        // the two middle samples to be symmetric.
 | 
						|
        hzv[mid] = fr * symsamples / 2.0;
 | 
						|
        hzv[mid - 1] = fr * symsamples / 2.0;
 | 
						|
    }
 | 
						|
 | 
						|
    // repeat first and last symbols
 | 
						|
    for (int i = 0; i < symsamples; i++)
 | 
						|
    {
 | 
						|
        hzv[i] = hzv[i + symsamples];
 | 
						|
        hzv[symsamples * (symbols.size() + 1) + i] = hzv[symsamples * symbols.size() + i];
 | 
						|
    }
 | 
						|
 | 
						|
    // run the per-sample frequency vector through
 | 
						|
    // the gaussian filter.
 | 
						|
    int half = gwin.size() / 2;
 | 
						|
    std::vector<float> o(hzv.size());
 | 
						|
    for (int i = 0; i < (int)o.size(); i++)
 | 
						|
    {
 | 
						|
        float sum = 0;
 | 
						|
        for (int j = 0; j < (int)gwin.size(); j++)
 | 
						|
        {
 | 
						|
            int k = i - half + j;
 | 
						|
            if (k >= 0 && k < (int)hzv.size())
 | 
						|
            {
 | 
						|
                sum += hzv[k] * gwin[j];
 | 
						|
            }
 | 
						|
        }
 | 
						|
        o[i] = sum;
 | 
						|
    }
 | 
						|
 | 
						|
    // drop repeated first and last symbols
 | 
						|
    std::vector<float> oo(symsamples * symbols.size());
 | 
						|
    for (int i = 0; i < (int)oo.size(); i++)
 | 
						|
    {
 | 
						|
        oo[i] = o[i + symsamples];
 | 
						|
    }
 | 
						|
 | 
						|
    // now oo[i] contains the frequency for the i'th sample.
 | 
						|
 | 
						|
    std::vector<float> v(symsamples * symbols.size());
 | 
						|
    float theta = phase0;
 | 
						|
    for (int i = 0; i < (int)v.size(); i++)
 | 
						|
    {
 | 
						|
        v[i] = cos(theta);
 | 
						|
        float hz = oo[i];
 | 
						|
        theta += 2 * M_PI / (rate / hz);
 | 
						|
    }
 | 
						|
 | 
						|
    return v;
 | 
						|
}
 | 
						|
 | 
						|
const std::string WHITESPACE = " \n\r\t\f\v";
 | 
						|
 | 
						|
std::string ltrim(const std::string &s)
 | 
						|
{
 | 
						|
    size_t start = s.find_first_not_of(WHITESPACE);
 | 
						|
    return (start == std::string::npos) ? "" : s.substr(start);
 | 
						|
}
 | 
						|
 | 
						|
std::string rtrim(const std::string &s)
 | 
						|
{
 | 
						|
    size_t end = s.find_last_not_of(WHITESPACE);
 | 
						|
    return (end == std::string::npos) ? "" : s.substr(0, end + 1);
 | 
						|
}
 | 
						|
 | 
						|
std::string trim(const std::string &s) {
 | 
						|
    return rtrim(ltrim(s));
 | 
						|
}
 | 
						|
 | 
						|
} // namespace FT8
 |