rearranging frequencies translation

pull/23/head
Xael South 2021-02-15 20:18:07 +00:00
rodzic f9f42caafb
commit 139633fa0b
2 zmienionych plików z 52 dodań i 61 usunięć

Wyświetl plik

@ -103,21 +103,24 @@ Run length algorithm works well with a few mode C1 devices I had around me, but
S1 mode datagrams can now be received! You have to start rtl_wmbus at 868.3MHz with
* rtl_sdr -f 868.3M -s 1600000 - 2>/dev/null | build/rtl_wmbus
time2 clock recovery method has to be activated in this case (it's active by default). It works very well because of Manchester coding.
Last but not least, you can try to receive all datagrams (S1, T1, C1) _simultaneously_:
* rtl_sdr -f 868.7M -s 1600000 - 2>/dev/null | build/rtl_wmbus -s
* rtl_sdr -f 868.625M -s 1600000 - 2>/dev/null | build/rtl_wmbus -s
Notice in the last line:
* "-s": which is needed to inform rtl_wmbus about required frequency translation
* "868.7M": the new frequency to receive at
* "868.625M": the new frequency to receive at
rtl_wmbus will then shift all frequencies
* by 250kHz to new center frequency at 868.95Mhz (T1 and C1)
* by 400kHz to new center frequency at 868.3Mhz (S1)
* by +325kHz to new center frequency at 868.95Mhz (T1 and C1)
* by -325kHz to new center frequency at 868.3Mhz (S1)
I have tested this so far and can confirm that it works for T1/C1 and S1. Thanks to alalons for providing me with bitstreams!
Optimization on frequencies translation by rearranging compute steps implemented as proposed by alalons.
Alalons (have I thanked you already?!) proposed an speed optimized arctan function. Performance gain is notable (factor ~2) but could reduce sensitivity slightly on receiving C1 mode datagrams. A speed optimized arctan version can be activated by "-a" in the program options.
License
-------

Wyświetl plik

@ -371,10 +371,10 @@ static inline float polar_discriminator_t1_c1(float i, float q)
static inline float polar_discriminator_t1_c1_inaccurate(float i, float q)
{
// We are using only complex part of the phase difference
// We are going to use only complex part of the phase difference
// so avoid unnecesary computation of real part. The math behind:
// cargf = atan (delta_phi_imag / delta_phi_real) / pi;
// In the formula only the sign is of interest.
// In the formula only the sign is of interest - we compute delta_phi_imag only.
static float i_last, q_last;
@ -407,10 +407,10 @@ static inline float polar_discriminator_s1(float i, float q)
static inline float polar_discriminator_s1_inaccurate(float i, float q)
{
// We are using only complex part of the phase difference
// We are going to use only complex part of the phase difference
// so avoid unnecesary computation of real part. The math behind:
// cargf = atan (delta_phi_imag / delta_phi_real) / pi;
// In the formula only the sign is of interest.
// In the formula only the sign is of interest - we compute delta_phi_imag only.
static float i_last, q_last;
const float delta_phi_imag = i_last*q - i*q_last;
@ -751,78 +751,68 @@ static void process_options(int argc, char *argv[])
}
}
static float complex *LUT_FREQUENCY_TRANSLATION_PLUS = NULL;
static float complex *LUT_FREQUENCY_TRANSLATION_MINUS = NULL;
#define FREQ_STEP_KHZ (50)
static float *LUT_FREQUENCY_TRANSLATION_PLUS_COSINE = NULL;
static float *LUT_FREQUENCY_TRANSLATION_PLUS_SINE = NULL;
#define FREQ_STEP_KHZ (25)
/* fs_kHz is the sample rate in kHz. */
void setup_lookup_tables_for_frequency_translation(int fs_kHz)
static void setup_lookup_tables_for_frequency_translation(int fs_kHz)
{
const int ft_kHz = FREQ_STEP_KHZ;
const size_t n_max = fs_kHz/ft_kHz;
free(LUT_FREQUENCY_TRANSLATION_PLUS);
LUT_FREQUENCY_TRANSLATION_PLUS = malloc(n_max *sizeof(LUT_FREQUENCY_TRANSLATION_PLUS[0]));
if (!LUT_FREQUENCY_TRANSLATION_PLUS) exit(EXIT_FAILURE);
free(LUT_FREQUENCY_TRANSLATION_PLUS_COSINE);
LUT_FREQUENCY_TRANSLATION_PLUS_COSINE = malloc(n_max *sizeof(LUT_FREQUENCY_TRANSLATION_PLUS_COSINE[0]));
if (!LUT_FREQUENCY_TRANSLATION_PLUS_COSINE) exit(EXIT_FAILURE);
free(LUT_FREQUENCY_TRANSLATION_MINUS);
LUT_FREQUENCY_TRANSLATION_MINUS = malloc(n_max *sizeof(LUT_FREQUENCY_TRANSLATION_MINUS[0]));
if (!LUT_FREQUENCY_TRANSLATION_MINUS) exit(EXIT_FAILURE);
free(LUT_FREQUENCY_TRANSLATION_PLUS_SINE);
LUT_FREQUENCY_TRANSLATION_PLUS_SINE = malloc(n_max *sizeof(LUT_FREQUENCY_TRANSLATION_PLUS_SINE[0]));
if (!LUT_FREQUENCY_TRANSLATION_PLUS_SINE) exit(EXIT_FAILURE);
for (size_t n = 0; n < n_max; n++)
{
const double phi = (2. * M_PI * (ft_kHz * n)) / fs_kHz;
const float a = cosf(phi);
const float b = sinf(phi);
LUT_FREQUENCY_TRANSLATION_PLUS[n] = a - b * _Complex_I;
LUT_FREQUENCY_TRANSLATION_MINUS[n] = a + b * _Complex_I;
LUT_FREQUENCY_TRANSLATION_PLUS_COSINE[n] = cosf(phi);
LUT_FREQUENCY_TRANSLATION_PLUS_SINE[n] = -sinf(phi); // Minus sinf!
}
}
/* Positive frequencies shift: ft = 250kHz the signal will shift from 868.7M right to 868.95M. */
static void shift_freq_250(float *i, float *q, int fs_kHz)
/* Positive frequencies shift: ft = +325kHz the signal will shift from 868.625M right to 868.95M.
Negative frequencies shift: ft = -325kHz the signal will shift from 868.625M right to 868.3M. */
static void shift_freq_plus_minus325(float *iplus, float *qplus, float *iminus, float *qminus, int fs_kHz)
{
const int ft = 250;
static size_t n = 0;
const size_t n_max = fs_kHz/FREQ_STEP_KHZ;
#if 0
const float complex freq_shift = cosf(2.*M_PI*(ft*n)/fs_kHz) - sin(2.*M_PI*(ft*n)/fs_kHz) * _Complex_I;
n++;
#else
const float complex freq_shift = LUT_FREQUENCY_TRANSLATION_PLUS[n];
n += ft/FREQ_STEP_KHZ;
#endif
//fprintf(stdout, "%ld, %ld: %f, %f\n", n, n_max, crealf(freq_shift), cimagf(freq_shift));
if (n >= n_max) n -= n_max;
float complex s = (*i) + (*q) * _Complex_I;
s = s * freq_shift;
*i = crealf(s);
*q = cimagf(s);
}
/* Negative frequencies shift: ft = 400kHz the signal will shift from 868.7M right to 868.3M. */
static void shift_freq_minus400(float *i, float *q, int fs_kHz)
{
const int ft = 400;
const int ft = 325;
static size_t n = 0;
const size_t n_max = fs_kHz/FREQ_STEP_KHZ;
#if 0
const float complex freq_shift = cosf(2.*M_PI*(ft*n)/fs_kHz) + sin(2.*M_PI*(ft*n)/fs_kHz) * _Complex_I;
n++;
#else
const float complex freq_shift = LUT_FREQUENCY_TRANSLATION_MINUS[n];
const float x = LUT_FREQUENCY_TRANSLATION_PLUS_COSINE[n];
const float z = LUT_FREQUENCY_TRANSLATION_PLUS_SINE[n];
n += ft/FREQ_STEP_KHZ;
#endif
//fprintf(stdout, "%ld, %ld: %f, %f\n", n, n_max, crealf(freq_shift), cimagf(freq_shift));
if (n >= n_max) n -= n_max;
float ix, iz, qx, qz;
float complex s = (*i) + (*q) * _Complex_I;
s = s * freq_shift;
// (i+Jq)*(x+Jz) =ix-qz + J(qx+iz) positive rotation
// (i+Jq)*(x-Jz) =ix+qz + J(qx-iz) negative rotation
// ix, qz, qx, iz are the same for boths shifts so we reuse them
// totaling 4 mul and 4 sums instead of 8 mul 4 sum.
// It works because iplus equals to iminus and qplus equals to qminus.
*i = crealf(s);
*q = cimagf(s);
ix = *iplus * x;
qx = *qplus * x;
iz = *iplus * z;
qz = *qplus * z;
// Symmetric positive shift.
*iplus = ix - qz;
*qplus = qx + iz;
// Symmetric negative shift .
*iminus = ix + qz;
*qminus = qx - iz;
}
int main(int argc, char *argv[])
@ -899,9 +889,7 @@ int main(int argc, char *argv[])
if (opts_s1_t1_c1_simultaneously)
{
shift_freq_250(&i_t1_c1_unfilt, &q_t1_c1_unfilt, fs_kHz);
shift_freq_minus400(&i_s1_unfilt, &q_s1_unfilt, fs_kHz);
shift_freq_plus_minus325(&i_t1_c1_unfilt, &q_t1_c1_unfilt, &i_s1_unfilt, &q_s1_unfilt, fs_kHz);
}
// Low-Pass-Filtering before decimation is necessary, to ensure
@ -1043,8 +1031,8 @@ int main(int argc, char *argv[])
}
if (input != stdin) fclose(input);
free(LUT_FREQUENCY_TRANSLATION_PLUS);
free(LUT_FREQUENCY_TRANSLATION_MINUS);
free(LUT_FREQUENCY_TRANSLATION_PLUS_COSINE);
free(LUT_FREQUENCY_TRANSLATION_PLUS_SINE);
return EXIT_SUCCESS;
}