Merge branch 'next' into next_memless

4 files changed, 101 insertions, 7 deletions

diff --git a/doc/example.ini b/doc/example.ini
index 0659813..f8cec36 100644
--- a/doc/example.ini
+++ b/doc/example.ini
@@ -141,13 +141,16 @@ output=uhd
 
 [fileoutput]
 ; Two output formats are supported: In the default mode,
-; the file output writes I/Q float values (i.e. complex float)
-; to the file. The I and Q samples can take values up to
-; 100000 in absolute magnitude with gainmode FIX.
-; With gainmode VAR, they should never exceed 50000.
-; With gainmode MAX, thet are limited to 32767.
+; the file output writes I/Q float values (i.e. complex
+; float) to the file. The I and Q samples can take values up
+; to 810000 in absolute magnitude with gainmode FIX. With
+; gainmode VAR and FIX, they should never exceed 50000.
 ;format=complexf
 ;
+; When the format is set to complexf_normalised the maximal
+; amplitude written to a file, is the digital gain +-10%.
+;format=complexf_normalised
+;
 ; When the format is set to s8, the output writes I/Q 8-bit
 ; signed integers, where the magnitude is multiplied by 128/50000
 ; effectively mapping the gainmode VAR range of -50000 -- 50000
diff --git a/doc/time-freq-plot.py b/doc/time-freq-plot.py
new file mode 100755
index 0000000..6ece564
--- /dev/null
+++ b/doc/time-freq-plot.py
@@ -0,0 +1,77 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+#
+# Print scope and spectrum from ODR-DabMod I/Q file
+#
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Matthias P. Braendli
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+#    The above copyright notice and this permission notice shall be included in
+#    all copies or substantial portions of the Software.
+#
+#    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+#    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+#    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+#    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+#    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+#    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+#    THE SOFTWARE.
+import sys
+import matplotlib.pyplot as plt
+import numpy as np
+
+rate=2048000
+# T = 1/2048000 s
+# NULL symbol is 2656 T (about 1.3ms) long.
+T_NULL = 2656
+# Full transmission frame in TM1 is 96ms = 196608 T.
+T_TF = 196608
+
+num_skip_samples = 8 * T_TF
+num_analyse_samples = 2 * T_TF
+
+if len(sys.argv) < 2:
+    print("Specify .iq file name")
+    print("Expected format: complex float I/Q, 2048000 Sps")
+    print("The input file must contain at least 10 transmission frames,")
+    print("i.e. {} samples = {} seconds".format(T_TF * 10, T_TF * 10.0 / rate))
+    sys.exit(1)
+
+fd = open(sys.argv[1], 'rb')
+
+# The IQ files potentially have zero samples in the beginning, we need
+# to skip a few transmission frames
+
+
+source_data = np.fromfile(file=fd, dtype=np.complex64, count=num_skip_samples + num_analyse_samples)
+
+print("Read in {} samples".format(len(source_data)))
+
+source_data = source_data[num_skip_samples:]
+source_data_time = np.linspace(0, num_analyse_samples/rate, len(source_data))
+
+print("Signal power: {} of {} samples".format(np.sum(np.abs(source_data**2)), len(source_data)))
+
+fft_size = 4096
+
+plt.figure(figsize=(10,8))
+plt.subplot(211)
+plt.title("Real part of signal")
+plt.plot(source_data_time, np.real(source_data))
+
+signal_spectrum = np.abs(np.fft.fftshift(np.fft.fft(source_data[T_NULL:], fft_size)))
+freqs = np.fft.fftshift(np.fft.fftfreq(fft_size, d=1./rate))
+
+plt.subplot(212)
+plt.title("Spectrum of {} samples after the NULL symbol".format(fft_size))
+plt.semilogy(freqs, signal_spectrum)
+
+plt.show()
diff --git a/src/DabMod.cpp b/src/DabMod.cpp
index 301b078..4e4cdab 100644
--- a/src/DabMod.cpp
+++ b/src/DabMod.cpp
@@ -81,6 +81,11 @@
  */
 static const float normalise_factor = 50000.0f;
 
+//Empirical normalisation factors used to normalise the samples to amplitude 1.
+static const float normalise_factor_file_fix = 81000.0f;
+static const float normalise_factor_file_var = 46000.0f;
+static const float normalise_factor_file_max = 46000.0f;
+
 typedef std::complex<float> complexf;
 
 using namespace std;
@@ -183,6 +188,15 @@ static shared_ptr<ModOutput> prepare_output(
         if (s.fileOutputFormat == "complexf") {
             output = make_shared<OutputFile>(s.outputName);
         }
+        if (s.fileOutputFormat == "complexf_normalised") {
+            if (s.gainMode == GainMode::GAIN_FIX)
+                s.normalise = 1.0 / normalise_factor_file_fix;
+            else if (s.gainMode == GainMode::GAIN_MAX)
+                s.normalise = 1.0 / normalise_factor_file_max;
+            else if (s.gainMode == GainMode::GAIN_VAR)
+                s.normalise = 1.0 / normalise_factor_file_var;
+            output = make_shared<OutputFile>(s.outputName);
+        }
         else if (s.fileOutputFormat == "s8") {
             // We must normalise the samples to the interval [-127.0; 127.0]
             s.normalise = 127.0f / normalise_factor;
diff --git a/src/Resampler.cpp b/src/Resampler.cpp
index ee2b865..8786e91 100644
--- a/src/Resampler.cpp
+++ b/src/Resampler.cpp
@@ -76,9 +76,9 @@ Resampler::Resampler(size_t inputRate, size_t outputRate, size_t resolution) :
     PDEBUG(" FFT size in: %zu, FFT size out: %zu\n", myFftSizeIn, myFftSizeOut);
 
     if (myFftSizeIn > myFftSizeOut) {
-        myFactor = 1.0f / myFftSizeIn;
+        myFactor = 1.0f / myFftSizeIn * outputRate / inputRate;
     } else {
-        myFactor = 1.0f / myFftSizeOut;
+        myFactor = 1.0f / myFftSizeOut * outputRate / inputRate;
     }
 
     myWindow = (float*)memalign(16, myFftSizeIn * sizeof(float));