Hałas kwantyzacji dla spójnego pobierania próbek - szum fazowy?

Aktualizacja: zobacz dodane przemyślenia na dole tego postu.

---

W ogólnych warunkach próbkowania, które nie są ograniczone tym, co opisano poniżej (sygnał niezwiązany z zegarem próbkowania), szum kwantyzacji jest często szacowany jako równomierny rozkład na jednym poziomie kwantyzacji. Kiedy dwa ADC są połączone ze ścieżkami I i Q w celu utworzenia próbkowania złożonego sygnału, szum kwantyzacji ma zarówno składową amplitudy, jak i szumu fazowego, jak symulowano poniżej. Jak pokazano, szum ten ma rozkład trójkątny, gdy składowe I i Q przyczyniają się w równym stopniu do amplitudy i fazy, na przykład gdy sygnał jest pod kątem 45 °, i jednolity, gdy sygnał jest na osi. Jest to oczekiwane, ponieważ szum kwantyzacji dla każdego I i Q jest nieskorelowany, więc rozkłady będą się zwoływać, gdy oba przyczynią się do wyniku wyjściowego.

Zadawane jest pytanie, czy ten rozkład szumu fazowego zmienia się znacząco w przypadkach spójnego próbkowania (zakładając, że sam zegar próbkowania ma szum fazowy, który jest znacznie lepszy, więc nie jest czynnikiem)? W szczególności staram się zrozumieć, czy spójne próbkowanie znacznie zmniejszy szum fazowy związany z kwantyzacją. Miałoby to bezpośrednie zastosowanie do generowania sygnału zegarowego, w którym spójność byłaby łatwa do utrzymania.

Rozważ zarówno sygnały rzeczywiste (jeden ADC), jak i sygnały złożone (dwa ADC; jeden dla I i jeden dla Q, razem opisując jedną złożoną próbkę). W przypadku sygnałów rzeczywistych wejściem jest fala sinusoidalna w pełnej skali, a składnik fazowy pochodzi z sygnału analitycznego; jitter związany ze zmianami w przejściu przez zero tonu sinusoidalnego byłby przykładem wynikowego szumu fazowego dla sygnału rzeczywistego. W przypadku sygnałów złożonych sygnał wejściowy jest w pełnej skali$Ae^{j \omega t}$, gdzie rzeczywistymi i urojonymi elementami byłyby fale sinusoidalne w pełnej skali.

Jest to związane z tym pytaniem, w którym spójne próbkowanie jest dobrze opisane, ale szum fazowy nie został konkretnie wspomniany:

Spójne pobieranie próbek i rozkład hałasu kwantyzacji

Aby jaśniej opisać indukowane składowe szumu AM i PM, dodałem poniższą grafikę dla przypadku złożonej kwantyzacji pokazującej złożony wektor w ciągłym czasie w danym momencie próbkowania i związaną z nim skwantowaną próbkę jako czerwoną kropkę, zakładając liniowy równomierny rozkład poziomów kwantyzacji rzeczywistych i urojonych części sygnału.

Zbliżenie na miejsce, w którym występuje kwantyzacja na powyższej grafice, aby zilustrować indukowany błąd amplitudy i błąd fazy:

W ten sposób otrzymano dowolny sygnał

$$\begin{align} s(t) &= a(t) e^{j\omega t} \\ &= a(t) \cos(\omega t) + j a(t) \sin(\omega t) \\ &= i(t) + j q(t) \\ \end{align}$$

Skwantowany sygnał jest najbliższym punktem odległości podanym przez

$$s_k = i_k+ j q_k$$

Gdzie $i_k$ i $q_k$ reprezentują skwantowane poziomy I i Q, każdy zamapowany zgodnie z:

$$\mathcal{Q}\{x\} = \Delta \Bigl \lfloor \frac{x}{\Delta}+\tfrac{1}{2} \Bigr \rfloor$$

Gdzie $\lfloor (\cdot) \rfloor$reprezentuje funkcję podłogi , oraz$\Delta$ reprezentuje dyskretny poziom kwantyzacji.

$$\begin{align} i_k = \mathcal{Q}\{i(t_k)\} \\ q_k = \mathcal{Q}\{q(t_k)\} \\ \end{align}$$

Błąd amplitudy wynosi $|s(t_k)|-|s_k|$ gdzie $t_k$ jest to czas $s(t)$ próbkowano w celu wygenerowania $s_k$.

Błąd fazy wynosi $\arg\{s(t_k)\} - \arg\{s_k\} = \arg\{s(t_k) \cdot (s_k)^*\}$ gdzie * reprezentuje koniugat złożony.

Pytanie do tego postu brzmi: jaka jest natura elementu fazy, gdy zegar próbkowania jest współmierny z (wielokrotnością całkowitą) sygnałem wejściowym?

Aby pomóc, oto niektóre symulowane rozkłady błędów amplitudy i fazy dla złożonego przypadku kwantyzacji z kwantyzacją 6 bitów na I i Q. Dla tych symulacji zakłada się, że rzeczywista „prawda” sygnału jest równie prawdopodobne w dowolnym miejscu kwantyzacji sektor zdefiniowany jako siatka pokazana na powyższym schemacie. Zwróć uwagę, gdy sygnał znajduje się wzdłuż jednej z ćwiartek (zarówno I, jak i Q), rozkład jest jednolity, jak oczekiwano w pojedynczym przypadku ADC z rzeczywistymi sygnałami. Ale gdy sygnał jest ustawiony pod kątem 45 °, rozkład jest trójkątny. Ma to sens, ponieważ w takich przypadkach sygnał ma równy udział I i Q, z których każdy jest nieskorelowanym rozkładem jednorodnym; więc te dwa rozkłady są trójwymiarowe.

Po obróceniu wektora sygnału do 0 ° histogramy wielkości i kąta są znacznie bardziej jednolite, zgodnie z oczekiwaniami:

---

Aktualizacja: Ponieważ nadal potrzebujemy odpowiedzi na konkretne pytanie (odpowiedź Olli poniżej zawiera dobre wyjaśnienie właściwości hałasu, który doprowadził do mojej aktualizacji gęstości trójkątnych i jednolitych hałasów, ale charakterystykę hałasu fazowego pod spójne warunki próbkowania są wciąż nieuchwytne), oferuję następujące przemyślenia, które mogą pobudzić rzeczywistą odpowiedź lub dalszy postęp (zauważ, że są to przemyślenia, które mogą być często mylone, ale w celu uzyskania odpowiedzi, której jeszcze nie mam):

Należy zauważyć, że w spójnych warunkach próbkowania częstotliwość próbkowania jest całkowitą wielokrotnością częstotliwości wejściowej (i również fazą zablokowaną). Oznacza to, że zawsze będzie liczba całkowita próbek, gdy raz obrócimy się przez płaszczyznę złożoną dla złożonego sygnału i próbkowania, lub liczba całkowita dla jednego cyklu sinusoidy dla rzeczywistego sygnału i próbkowania (pojedynczy ADC).

I jak opisano, zakładamy, że sam zegar próbkowania jest znacznie lepszy, więc nie jest uważany za wkład. Dlatego próbki za każdym razem lądują w dokładnie tym samym miejscu.

Biorąc pod uwagę przypadek rzeczywistego sygnału, gdybyśmy zajmowali się tylko przejściem przez zero przy określaniu szumu fazowego, wynikiem spójnego próbkowania byłoby jedynie ustalone, ale spójne przesunięcie opóźnienia (chociaż zbocza narastające i opadające mogą mieć różne opóźnienia gdy spójność jest nieparzystą liczbą całkowitą). Najwyraźniej w złożonym przypadku próbkowania mamy do czynienia z szumem fazowym przy każdej próbce i podejrzewam, że byłby to również ten sam przypadek rzeczywisty (podejrzewam, że opóźnienie próbki w dowolnym momencie od „prawdy” byłoby składnik szumu fazowego, ale potem się mylę, jeśli podwójnie liczę, jaka jest również różnica amplitudy ...) Jeśli mam czas, zasymuluję to, ponieważ wszystkie zniekształcenia pojawią się przy całkowitych harmonicznych sygnału wejściowego, biorąc pod uwagę powtarzający się wzór na jeden cykl, a test fazy w stosunku do amplitudy byłby względną fazą harmonicznych w stosunku do fundamentalnej - to, co byłoby interesujące zobaczyć za pomocą symulacji lub obliczeń, to czy te harmoniczne (które dla prawdziwego sygnału miałyby złożone sprzężone odpowiedniki) sumują się w kwadraturze z podstawową lub w fazie, a zatem wykazano, że jest to cały szum fazowy, cały hałas amplitudowy lub złożony z nich obu. (Różnica między parzystą liczbą próbek a liczbą nieparzystą może prawdopodobnie mieć na to wpływ).

W przypadku kompleksu grafika Olli'ego, która została wykonana z proporcjonalną liczbą próbek, może dodać więcej wglądu, jeśli pokazał lokalizację próbki na „prawdzie” związanej z każdą pokazaną skwantowaną próbką. Znów widzę możliwość interesującej różnicy, jeśli istnieje nieparzysta lub parzysta liczba próbek (jego grafika była parzysta i obserwuję wynikającą z niej symetrię, ale nie widzę dalej tego, co może zrobić z szumem fazowym w stosunku do amplitudy). Wydaje mi się jednak jasne, że składowe szumu zarówno w rzeczywistych, jak i złożonych przypadkach będą istniały tylko przy całkowitych harmonicznych częstotliwości podstawowej, gdy próbkowanie jest spójne. Więc nawet jeśli szum fazowy może nadal istnieć, tak jak podejrzewam, to jego lokalizacja przy harmonicznych liczbach całkowitych jest o wiele bardziej sprzyjająca eliminacji przez kolejne filtrowanie.

(Uwaga: ma to zastosowanie do generowania referencyjnych sygnałów zegarowych o wysokiej czystości widmowej).

Mam wątpliwości dotyczące (Edytuj: zostało to później usunięte z pytania):

Można rozsądnie założyć, że rozkład tych składowych szumu AM i PM jest jednolity, o ile sygnał wejściowy nie jest skorelowany z zegarem próbkowania

Rozważ sygnał:

$$\operatorname{signal}(t) = \cos(t) + j\sin(t)$$ i jego kwantyzacja: $$\operatorname{quantized\_signal}(t) = \frac{\operatorname{round}\big(N\cos(t)\big)}{N} + j\times\frac{\operatorname{round}\big(N\sin(t)\big)}{N}$$

dla kroku kwantyzacji wynoszącego $1/N$ obu składników I i Q (masz $N = 5$ na twojej figurze).

Ryc. 1. Ślad sygnału (niebieska linia) i jego kwantyzacja (czarne kropki) oraz morfing między nimi, aby zobaczyć, w jaki sposób kwantyzowane są różne części sygnału, dla $N=5$. „Morphing” to po prostu zestaw dodatkowych wykresów parametrycznych$a\operatorname{signal}(t) + (1-a)\operatorname{quantized\_signal}(t)$ w $a = \left[\frac{1}{5}, \frac{2}{5},\frac{3}{5},\frac{4}{5}\right].$

Błąd w fazie spowodowany błędem kwantyzacji wynosi:

$$\operatorname{phase\_error}(t) = \operatorname{atan}\Big(\operatorname{Im}\big(\operatorname{quantized\_signal}(t)\big), \operatorname{Re}\big(\operatorname{quantized\_signal}(t)\big)\Big)\\- \operatorname{atan}\Big(\operatorname{Im}\big(\operatorname{signal}(t)\big), \operatorname{Re}\big(\operatorname{signal}(t)\big)\Big) \\= \operatorname{atan}\Big(\operatorname{round}\big(N\sin(t)\big), \operatorname{round}\big(N\cos(t)\big)\Big) - \operatorname{atan}\big(N\sin(t), N\cos(t)\big) \\= \operatorname{atan}\Big(\operatorname{round}\big(N\sin(t)\big), \operatorname{round}\big(N\cos(t)\big)\Big) - \operatorname{mod}(t-\pi, 2\pi) + \pi$$

Odejmowanie owiniętych faz jest ryzykowne, ale w tym przypadku działa.

Rysunek 2. $\operatorname{phase\_error}(t)$ dla $N = 5$.

To elementarna funkcja liniowa. Wszystkie segmenty linii przekraczają poziom zerowy, ale kończą się na różnych innych poziomach. To znaczy, biorąc pod uwagę$t$ jako jednolita zmienna losowa, która w funkcji gęstości prawdopodobieństwa $\operatorname{phase\_error}(t),$wartości bliskie zeru są nadmiernie reprezentowane. Więc$\operatorname{phase\_error}(t)$ nie może mieć jednolitego rozkładu.

Biorąc pod uwagę rzeczywiste pytanie, patrząc na ryc. 1, wystarczająco wysoko $N$i na takiej częstotliwości złożonego sinusoidy, że podczas każdego okresu próbkowania sygnał obraca się poza kilkoma granicami kwantyzacji, błędy kwantyzacji w próbkach są w rzeczywistości ustaloną sekwencją liczb pseudolosowych pochodzących z dziwactw teorii liczb. Błędy zależą od częstotliwości i od$N,$a także w fazie początkowej, jeśli częstotliwość jest podwielokrotnością wielokrotności częstotliwości próbkowania, w którym to przypadku błąd kwantyzacji jest powtarzającą się sekwencją, która nie zawiera wszystkich możliwych wartości błędu kwantyzacji. W granicy dużych$N$rozkłady błędów I i Q są jednolite, a błędy fazy i wielkości są liczbami pseudolosowymi pochodzącymi z rozkładów zależnych od fazy sygnału. Istnieje zależność od fazy, ponieważ prostokątna siatka kwantyzacji ma orientację.

W granicy dużych $N,$błąd fazowy i błąd wielkości są prostopadłymi składowymi błędu zespolonego. Błąd wielkości można wyrazić proporcjonalnie do nieskończenie małego kroku kwantyzacji, a błąd fazy można wyrazić proporcjonalnie do$\arcsin$kroku kwantyzacji. W fazie sygnału$\alpha$ błąd wielkości jest w kierunku kątowym $\alpha$ i błąd fazy jest w kierunku kątowym $\alpha + \pi/2$. Złożony błąd kwantyzacji rozkłada się równomiernie w kwadracie kroku kwantyzacji zorientowanym wzdłuż osi I i Q, z narożami o współrzędnych wyrażonych proporcjonalnie do kroku kwantyzacji:

$$\big[(1/2, 1/2),\quad(-1/2, 1/2),\quad(-1/2, -1/2),\quad(1/2, -1/2)\big]$$

Obrót tych współrzędnych lub ich równoważne rzutowanie na proporcjonalny błąd fazy i proporcjonalny błąd błędu wielkości daje zarówno tę samą płaską górną, jak i kawałkową funkcję liniowej gęstości prawdopodobieństwa z węzłami:

$$\left[\frac{\cos(\alpha)}{2} - \frac{\sin(\alpha)}{2},\quad \frac{\cos(\alpha)}{2} + \frac{\sin(\alpha)}{2},\quad -\frac{\cos(\alpha)}{2} + \frac{\sin(\alpha)}{2},\quad -\frac{\cos(\alpha)}{2} - \frac{\sin(\alpha)}{2}\right] = \left[\sqrt{2}\cos(\alpha + \pi/4),\quad \sqrt{2}\sin(\alpha + \pi/4),\quad -\sqrt{2}\cos(\alpha + \pi/4),\quad -\sqrt{2}\sin(\alpha + \pi/4)\right]$$

Rysunek 3. Węzły współdzielonej częściowej liniowej płaskiej góry gęstości gęstości podatności (PDF) proporcjonalnego błędu fazy i błędu proporcjonalnej wielkości, biorąc pod uwagę kąt sygnału $\alpha$. W$\alpha \in \{-\pi, -\pi/2, 0, \pi/2, \pi\}$PDF jest prostokątny. Niektóre węzły łączą się również w$\alpha \in \{-3\pi/4, -\pi/4, \pi/4, 3\pi/4\}$ dając trójkątny plik PDF z najgorszym przypadkiem$N$ asymptotyczne oszacowanie 1) maksymalnego błędu bezwzględnej wielkości wynoszącego $\sqrt{2}/2$ kroki kwantyzacji i 2) maksymalny bezwzględny błąd fazy $\sqrt{2}/2$ czasy $\arcsin$ kroku kwantyzacji.

W fazach pośrednich plik PDF wygląda na przykład tak:

Rysunek 4. Udostępniony plik PDF pod adresem $\alpha = \pi/8.$

Jak sugeruje Dan, plik PDF jest również zbiorem prostokątnych plików PDF błędów I i Q rzutowanych na osie błędu wielkości i fazy. Szerokość jednego z rzutowanych plików PDF wynosi$|\cos(\alpha)|$, a szerokość drugiego wynosi $|\sin(\alpha)|$. Ich łączna wariancja to$\cos^2(\alpha)/12 + \sin^2(\alpha)/12 = 1/12,$ mundur ponad $\alpha$.

Mogą istnieć pewne „pseudolucky” kombinacje fazy początkowej i stosunek liczby wymiernej częstotliwości złożonej sinusoidy i częstotliwości próbkowania, które dają tylko mały błąd dla wszystkich próbek w powtarzającej się sekwencji. Ze względu na symetrię błędów widocznych na ryc. 1, w sensie maksymalnego błędu bezwzględnego częstotliwości te są korzystne, dla których liczba punktów odwiedzonych na okręgu jest wielokrotnością 2, ponieważ szczęście (niski błąd) jest potrzebne przy tylko połowa punktów. Błąd w pozostałych punktach jest duplikatem tego, co jest na pierwszych, z odwróceniem znaku. Co najmniej wielokrotności 6, 4 i 12 mają jeszcze większą przewagę. Nie jestem pewien, jaka jest dokładna zasada, ponieważ nie wydaje się, że chodzi o bycie wielokrotnością czegoś. To' chodzi o symetrie siatki połączone z arytmetyką modulo. Niemniej jednak błędy pseudolosowe są deterministyczne, więc wyczerpujące poszukiwanie ujawnia najlepsze aranżacje. Znalezienie najlepszych rozwiązań w sensie błędu bezwzględnego pierwiastkowego (RMS) jest najłatwiejsze:

Rysunek 5. Góra) Najniższe możliwe błędy bezwzględnej kwantyzacji RMS w złożonym oscylatorze IQ dla różnych głębokości bitów oscylatora, z wykorzystaniem kwadratowej siatki kwantyzacji . Kod źródłowy do wyczerpującego poszukiwania pseudolubnych rozwiązań znajduje się na końcu odpowiedzi. Dół) Detal, pokazujący do porównania (jasnoniebieski)$N\to\infty$ asymptotyczne oszacowanie bezwzględnego błędu kwantyzacji RMS, $\sqrt{1/6}/N,$ dla $N=2^k-1,$ gdzie $k+1$ to liczba bitów oscylatora.

Amplituda najbardziej widocznej częstotliwości błędu nigdy nie jest większa niż błąd bezwzględny RMS. Szczególnie dobrym wyborem dla 8-bitowego oscylatora$12$ punkty znajdujące się w przybliżeniu na okręgu jednostki:

$$\frac{\{(0, \pm112),\quad (\pm112, 0),\quad (\pm97, \pm56),\quad (\pm56, \pm97)\}}{112.00297611139371}$$

Dyskretna złożona sinusoida, która przechodzi przez te punkty na płaszczyźnie złożonej w rosnącym porządku kątowym, ma tylko zniekształcenie 5. harmoniczne, a przy $-91.5$ dB w porównaniu do podstawowego, co potwierdza kod źródłowy Octave na końcu odpowiedzi.

Aby uzyskać bezwzględny błąd kwantyzacji niskiego RMS, częstotliwości nie muszą przechodzić przez punkty w kolejności jak w przybliżonych fazach $[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]\cdot 2\pi/12$ na częstotliwość $1/12$razy częstotliwość próbkowania. Na przykład częstotliwość$5/12$ razy częstotliwość próbkowania przejdzie przez te same punkty, ale w innej kolejności: $[0, 5, 10, 3, 8, 1, 6, 11, 4, 9, 2, 7]\cdot 2\pi/12$. Myślę, że to działa tak, jak działa, ponieważ 5 i 12 są chronione prawem autorskim .

Jeśli chodzi o możliwe idealne ustawienia, błąd może wynosić dokładnie zero we wszystkich punktach, jeśli częstotliwość sinusoidy stanowi jedną czwartą częstotliwości próbkowania (przyrost fazy $\pi/2$na próbkę). Na kwadratowej siatce nie ma innych tak doskonałych aranżacji . Na sześciokątnej siatce lub na kwadratowej prostokątnej siatce z jedną z osi I lub Q rozciągniętymi współczynnikiem$\sqrt{3}$ (przy czym odpowiada to co drugiemu rzędowi na siatce o strukturze plastra miodu), przyrost fazy o $\pi/3$na próbkę działałoby idealnie. Takie skalowanie można wykonać w domenie analogowej. Zwiększa to liczbę osi symetrii siatki, co powoduje w większości korzystne zmiany w pseudolubnych układach:

Rysunek 6. Najniższe możliwe błędy bezwzględnej kwantyzacji RMS w złożonym oscylatorze IQ dla różnych głębokości bitów oscylatora, przy użyciu prostokątnej siatki kwantyzacji z jedną osią skalowaną przez$\sqrt{3}$.

Warto zauważyć, że w przypadku 8-bitowego oscylatora z 30 punktami na okręgu najmniejszy możliwy błąd bezwzględny RMS wynosi -51,3 dB na kwadratowej siatce i -62,5 dB na kwadratowej prostokątnej siatce, gdzie najniższy błąd bezwzględny RMS błąd pseudolucky ma błąd:

Rysunek 7. Wartości błędu na płaszczyźnie IQ 8-bitowej pseudolubnej sekwencji długości 30 wykorzystują osie symetrii znalezione w siatce kwantyzacji rozciągniętej o współczynnik $\sqrt{3}$poziomo. Punkty pochodzą z zaledwie trzech pseudolubnych liczb zespolonych, odwróconych wokół osi symetrii.

Nie mam praktycznego doświadczenia z sygnałami zegara IQ, więc nie jestem pewien, co się liczy. Przy generowaniu sygnału zegarowego za pomocą przetwornika cyfrowo-analogowego (DAC) podejrzewam, że o ile nie zostaną zastosowane dobre układy pseudolucky, lepiej jest mieć niższy poziom białego szumu niż harmoniczne widmo hałasu o wyższym impulsy pochodzące z powtarzającej się sekwencji błędu kwantyzacji (patrz Spójne próbkowanie i rozkład hałasu kwantyzacji ). Te widmowe skoki, podobnie jak biały szum, mogą przeciekać przez pasożytniczą pojemność i mieć niepożądane skutki w innych częściach systemu lub wpływać na kompatybilność elektromagnetyczną (EMC) urządzenia. Analogicznie, technologia rozproszonego widma poprawia EMC poprzez zamianę widmowych pików na niższy poziom szumu.

Poniżej znajduje się kod źródłowy do wyczerpującego wyszukiwania pseudolucky aranżacji w C ++. Możesz uruchomić go na noc, aby znaleźć najlepsze ustawienia dla co najmniej 16-bitowych oscylatorów$1 \le M \le 100$.

// Compile with g++ -O3 -std-c++11

#include <stdio.h>
#include <math.h>
#include <complex>
#include <float.h>
#include <algorithm>

// N = circle size in quantization steps
const int maxN = 127;
// M = number of points on the circle
const int minM = 1; 
const int maxM = 100;
const int stepM = 1;
// k = floor(log2(N))
const int mink = 2;
const double IScale = 1; // 1 or larger please, sqrt(3) is very lucky, and 1 means a square grid

typedef std::complex<double> cplx;

struct Arrangement {
  int initialI;
  int initialQ;
  cplx fundamentalIQ;
  double fundamentalIQNorm;
  double cost;
};

int main() {
  cplx rotation[maxM+1];
  cplx fourierCoef[maxM+1];
  double invSlope[maxM+1];
  Arrangement bestArrangements[(maxM+1)*(int)(floor(log2(maxN))+1)];
  const double maxk(floor(log2(maxN)));
  const double IScaleInv = 1/IScale;
  for (int M = minM; M <= maxM; M++) {
    rotation[M] = cplx(cos(2*M_PI/M), sin(2*M_PI/M));
    invSlope[M] = tan(M_PI/2 - 2*M_PI/M)*IScaleInv;
    for (int k = 0; k <= maxk; k++) {
      bestArrangements[M+(maxM+1)*k].cost = DBL_MAX;
      bestArrangements[M+(maxM+1)*k].fundamentalIQNorm = 1;
    }
  }
  for (int M = minM; M <= maxM; M += stepM) {
    for (int m = 0; m < M; m++) {
      fourierCoef[m] = cplx(cos(2*M_PI*m/M), -sin(2*M_PI*m/M))/(double)M;
    }
    for (int initialQ = 0; initialQ <= maxN; initialQ++) {
      int initialI(IScale == 1? initialQ : 0);
      initialI = std::max(initialI, (int)floor(invSlope[M]*initialQ));
      if (initialQ == 0 && initialI == 0) {
    initialI = 1;
      }
      for (; initialI*(int_least64_t)initialI  <= (2*maxN + 1)*(int_least64_t)(2*maxN + 1)/4 - initialQ*(int_least64_t)initialQ; initialI++) {
    cplx IQ(initialI*IScale, initialQ);
    cplx roundedIQ(round(real(IQ)*IScaleInv)*IScale, round(imag(IQ)));
        cplx fundamentalIQ(roundedIQ*fourierCoef[0].real());
    for (int m = 1; m < M; m++) {
      IQ *= rotation[M];
      roundedIQ = cplx(round(real(IQ)*IScaleInv)*IScale, round(imag(IQ)));
          fundamentalIQ += roundedIQ*fourierCoef[m];
    }
    IQ = fundamentalIQ;
    roundedIQ = cplx(round(real(IQ)*IScaleInv)*IScale, round(imag(IQ)));
    double cost = norm(roundedIQ-IQ);
    for (int m = 1; m < M; m++) {
      IQ *= rotation[M];
      roundedIQ = cplx(round(real(IQ)*IScaleInv)*IScale, round(imag(IQ)));
      cost += norm(roundedIQ-IQ);
    }
    double fundamentalIQNorm = norm(fundamentalIQ);
    int k = std::max(floor(log2(initialI)), floor(log2(initialQ)));
    //  printf("(%d,%d)",k,initialI);
    if (cost*bestArrangements[M+(maxM+1)*k].fundamentalIQNorm < bestArrangements[M+(maxM+1)*k].cost*fundamentalIQNorm) {
      bestArrangements[M+(maxM+1)*k] = {initialI, initialQ, fundamentalIQ, fundamentalIQNorm, cost};
    }
      }
    }
  }
  printf("N");
  for (int k = mink; k <= maxk; k++) {
    printf(",%d-bit", k+2);
  }
  printf("\n");
  for (int M = minM; M <= maxM; M += stepM) {
    printf("%d", M);
    for (int k = mink; k <= maxk; k++) {
      printf(",%.13f", sqrt(bestArrangements[M+(maxM+1)*k].cost/bestArrangements[M+(maxM+1)*k].fundamentalIQNorm/M));
    }
    printf("\n");
  }

  printf("bits,M,N,fundamentalI,fundamentalQ,I,Q,rms\n");
  for (int M = minM; M <= maxM; M += stepM) {
    for (int k = mink; k <= maxk; k++) {
      printf("%d,%d,%.13f,%.13f,%.13f,%d,%d,%.13f\n", k+2, M, sqrt(bestArrangements[M+(maxM+1)*k].fundamentalIQNorm), real(bestArrangements[M+(maxM+1)*k].fundamentalIQ), imag(bestArrangements[M+(maxM+1)*k].fundamentalIQ), bestArrangements[M+(maxM+1)*k].initialI, bestArrangements[M+(maxM+1)*k].initialQ, sqrt(bestArrangements[M+(maxM+1)*k].cost/bestArrangements[M+(maxM+1)*k].fundamentalIQNorm/M));
    }
  }
}

Przykładowe dane wyjściowe opisujące pierwszą przykładową sekwencję znalezioną za pomocą IScale = 1:

bits,M,N,fundamentalI,fundamentalQ,I,Q,rms
8,12,112.0029761113937,112.0029761113937,0.0000000000000,112,0,0.0000265717171

Przykładowe dane wyjściowe opisujące drugą przykładową sekwencję znalezioną za pomocą IScale = sqrt(3):

8,30,200.2597744568315,199.1627304588310,20.9328464782995,115,21,0.0007529202390

Kod oktawowy do testowania pierwszej przykładowej sekwencji:

x = [112+0i, 97+56i, 56+97i, 0+112i, -56+97i, -97+56i, -112+0i, -97-56i, -56-97i, 0-112i, 56-97i, 97-56i];
abs(fft(x))
20*log10(abs(fft(x)(6)))-20*log10(abs(fft(x)(2)))

Kod oktawowy do testowania drugiej przykładowej sekwencji:

x = exp(2*pi*i*(0:29)/30)*(199.1627304588310+20.9328464782995i);
y = real(x)/sqrt(3)+imag(x)*i;
z = (round(real(y))*sqrt(3)+round(imag(y))*i)/200.2597744568315;
#Error on IQ plane
star = z-exp(2*pi*i*(0:29)/30)*(199.1627304588310+20.9328464782995i)/200.2597744568315;
scatter(real(star), imag(star));
#Magnitude of discrete Fourier transform
scatter((0:length(z)-1)*2*pi/30, 20*log10(abs(fft(z))/abs(fft(z)(2)))); ylim([-120, 0]);
#RMS error:
10*log10((sum(fft(z).*conj(fft(z)))-(fft(z)(2).*conj(fft(z)(2))))/(fft(z)(2).*conj(fft(z)(2))))