Author Topic: Idea for improving LTDZ Spectrum Analyzer & Tracking Generator 35MHz-4.4GHz  (Read 24172 times)

0 Members and 1 Guest are viewing this topic.

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
Edit: Regarding 80kHz: Yet another point to consider is the frequency deviation between the PLL xtal and the MCU xtal. If the (normalized) deviation between the two oscillators is (say) 20ppm, this leads to a frequency error of 0.3Hz when the IF frequency is 15kHz, but at 80kHz IF, the error is already 1.6Hz. However, the filter bandwidth is independent of the IF frequency. This means that the frequency error to bandwidth ratio is higher at 80kHz than at 15kHz. Still no problem yet and 1.6Hz error is still negligible (at least for magnitude estimation). But when we reduce the filter bandwidth more and more (by vector accumulation of more and more samples), in order to get more and more dynamic range, then at some point the frequency error will begin to matter. I'm also unsure whether the frequency deviation of all LTDZ devices is as low as on yours.

Changing the LO offset is quite easy in the firmware. All that is required is changing this line in the source code:
Code: [Select]
/** Default RX LO frequency offset in Hz in NA mode when using new ADC dB mode */
#define CONFIG_RX_DEFAULT_NA_MODE_FREQ_OFFSET_Hz    (10*8000L)
and regenerating the sin/cos-tables for the given LO offset frequency with the windowing applied. For this I have GNU Octave/Matlab script.
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
In this example, the vector sums for I and Q can continue accumulating

This implies eventually that you implicitly repeat the shorter (975 samples) window function N times, in order to to form the total window function which spans all N*975 samples. And this N-fold repetition leads to the weird frequency response then. But although the frequency response of this filter looks weird, it can still work for the given purpose, if the signal spectrum contains only noise at the position of the high side lobe (which is likely granted here).

For example I do not have strong digital signal processing background, so it has been real pleasure to have you here explaining things from the theoretical and practical point of view.  :-+

I kind of understand what is going on here, and why my intuition may be wrong, as you explained above. What we are basically doing is performing one tap DFT for the 975 samples [using a quadrature correlator]. Since DFT assumes that the signal to be analyzed is repetitive, proper windowing needs to be applied (or the buffer length needs to be selected appropriately so that the buffer length is exact multiple of the DFT bin frequency).

On the other hand, if we repeat this same process N times for different signal sample sets (while maintaining the phase, which is true in this case), this process will become as coherent averaging (Lyons, Understanding Digital Signal Processing, 3rd. ed. Chapters 11.1 Coherent Averaging and 11.3 Averaging Multiple Fast Fourier Transforms). As I understand, this  process is similar to as if were are performing [synchronous] FFT averaging over multiple FFT transforms with zero time-domain overlap. At the end of the chapter Lyons states that " ... it [coherent FFT integration] only works for periodic time-domain signal sequences that have been obtained through careful synchronous sampling. Coherent integration of multiple FFT results is of no value in reducing spectral measurement noise for non-periodic real-world, information carrying signals". Since we are satisfying this property (we are sampling stationary sinusoidal wave from LO mixer output with synchronous sampling*), it is possible to use coherent DFT integration and obtain increase in SNR. Please correct me if I have understood this wrong.

Edit: OTOH, when I wrote "we are sampling stationary sinusoidal wave from LO mixer output with synchronous sampling", this may not be 100% true because there is a small frequency difference between the ADC sampling clock and LO mixer output frequency, although it will  remain the same through out the sampling process in practical terms. As there is a small frequency offset present, the phase will change a little between sample sets. So, I am not quite sure how much this will impact the SNR gain, as we are using qudrature correlator.
« Last Edit: June 17, 2021, 07:47:28 am by Kalvin »
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
About synchronous sampling: Would it make any sense to use an ADC sample buffer of, let's say, 975 samples for 80kHz LO offset, and average N*975 samples into another "working" buffer, and only after N averaged sample sets compute the DFT? If the phase stays the same in practical terms between sample sets, this would reduce the real-time processing requirements, and allow longer sampling times without need for longer sample buffers. If the phase changes too much between sample sets, the SNR gain will be reduced (because the sampling process is no longer synchronous).
 

Offline gf

  • Super Contributor
  • ***
  • Posts: 1114
  • Country: de
For example I do not have strong digital signal processing background

Neither do I -- just a little bit.

Quote
Since DFT assumes that the signal to be analyzed is repetitive, proper windowing needs to be applied (or the buffer length needs to be selected appropriately so that the buffer length is exact multiple of the DFT bin frequency).

You already "apply windowing" when you pick out a finite sub-sequence of consecutive samples from the infinite stream of samples that represent the discrete signal from t=-inf...inf.
The window function just defines the relative weighting of these samples then (rectangular = equally weighted).
Stricly speaking, "no window" does actually not exist, but there is always some window involved when you truncate an infinite signal to a finite number of samples.

Quote
On the other hand, if we repeat this same process N times for different signal sample sets (while maintaining the phase, which is true in this case), this process will become as coherent averaging (Lyons, Understanding Digital Signal Processing, 3rd. ed. Chapters 11.1 Coherent Averaging and 11.3 Averaging Multiple Fast Fourier Transforms). As I understand, this  process is similar to as if were are performing [synchronous] FFT averaging over multiple FFT transforms with zero time-domain overlap. At the end of the chapter Lyons states that " ... it [coherent FFT integration] only works for periodic time-domain signal sequences that have been obtained through careful synchronous sampling. Coherent integration of multiple FFT results is of no value in reducing spectral measurement noise for non-periodic real-world, information carrying signals". Since we are satisfying this property (we are sampling stationary sinusoidal wave from LO mixer output with synchronous sampling*), it is possible to use coherent DFT integration and obtain increase in SNR. Please correct me if I have understood this wrong.

Yes, you basically integrate over a the full window length, i.e. over a contiguous sequence of equally spaced samples. This implies that the ADC must not be stopped and restarted asynchronously for each chunk if this sequence of samples is acquired in multiple chunks. If you skip M samples somewhere in the midde, and insert M zeros instead, then coherency is not lost, but it has the same effect then as if the window function were zero at these samples. But the frequency response of such a window function may become weird, as you have seen.

Quote
Edit: OTOH, when I wrote "we are sampling stationary sinusoidal wave from LO mixer output with synchronous sampling", this may not be 100% true because there is a small frequency difference between the ADC sampling clock and LO mixer output frequency, although it will  remain the same through out the sampling process in practical terms. As there is a small frequency offset present, the phase will change a little between sample sets. So, I am not quite sure how much this will impact the SNR gain, as we are using qudrature correlator.

For the amplidue estimation it means that the IF frequency does not hit exactly the center of the bandpass filter response, therefore the amplitude estimate will be a little bit lower.

If you take two readings at times t, and t+dt, and want to measure their phase difference, then the frequency error integrates over the interval dt to a phase error. If the error is larger than negligible, then the frequency offset needs to be estimated in order to compensate the phase error. Btw, the two (independent) readings at times t and t+dt need to be captured coherently, too, i.e. TG keeps running, LO keeps running, and ADC keeps running as well, only the receiver input is switched between the two readings. Otherwise phase information were lost in between. Still the interval dt should be kept short, since the frequency error can only be estimated with limited accuracy.

Quote
About synchronous sampling: Would it make any sense to use an ADC sample buffer of, let's say, 975 samples for 80kHz LO offset, and average N*975 samples into another "working" buffer, and only after N averaged sample sets compute the DFT? If the phase stays the same in practical terms between sample sets, this would reduce the real-time processing requirements, and allow longer sampling times without need for longer sample buffers. If the phase changes too much between sample sets, the SNR gain will be reduced (because the sampling process is no longer synchronous).

The following two expressions are mathematically equivalent:
Code: [Select]
sum(sum(x .* repmat(lo,N,1) .* repmat(win,N,1)))
sum(sum(x) .* lo .* win)
[ where x is a Nx975 matrix, and win and lo are 1x975 row vectors ]

So I'd say this is at the end equivalent to your previous proposal, which was in turn equivalent to an N-fold repetition of the 975-tap windows function (which lead to the weird frequency response).
« Last Edit: June 17, 2021, 02:14:39 pm by gf »
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
Let's assume that we are sampling N samples into a buffer. We have selected N in such way that ADC sampling frequency / LO frequency over N will become integer in all practical terms (for Fs=72Mhz/6/14=857142.8571428572_Hz and  LO=80kHz: Fs/LO=10.71428571428, we could select N = 93 * 10.71428571428 = 996.42857142857 = 996, for example). Now, we will compute the DFT for the given LO frequency of 80kHz (using quadrature correlator with pre-computed sin/cos tables with Hann-window). The RMS for this DFT = sqrt(sum(Q)^2 + sum(I)^2). No weird spectrum here.

Now, let's extend this a little, and that we are sampling K*N consecutive samples into a buffer. We will then divide this buffer into K blocks, each size of N samples. Then, we will compute DFT and RMS for each individual blocks of N samples. Finally, we will combine the results from all K blocks into one RMS value. There should not be no weird spectrum here, because we are just averaging K * individual DFT results. Right?
« Last Edit: June 19, 2021, 09:41:56 am by Kalvin »
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
Let's continue now from my previous post.

We have chosen N is such a clever way that the buffer of N samples will always contain integer number of LO cycles, thus the sin/cos-lookup tables will also always contain integer number of sin/cos cycles. We have also selected N so that it is even. This means that there won't be any sin/cos nor LO phase discontinuation between the K blocks, thus we are essentially performing coherent sampling for the K buffers (ie. K*N samples).

We can now do the same thing as previously "We will then divide this buffer into K blocks, each size of N samples. Then, we will compute DFT and RMS for each individual blocks of N samples. Finally, we will combine the results from all K blocks into one RMS value. There should not be no weird spectrum here, because we are just averaging K * individual DFT results."

But, we can do better than that. Since we sampling in coherent way, we do not have to compute the RMS sum for each individual block of N samples, but we can integrate the RMS sum over all K blocks, which will improve the SNR. However, this should not produce weird spectrum either.

The idea here is that we need to have only one set of sin/cos tables of N samples (with Hann-windowing), and we can select the K to be from 1 to any integer number we want to have. If we want fast sweep rate, we just select K=1. If we want better dynamic range, will select K to be 2, 3, 4 ... which will result slower sweep rate but improve the dynamic range.

Alternatively, because we are using coherent sampling, we could have a buffer of N samples in which we maintain the sum for each N samples. After be have captured all K blocks (K*N samples in total), we could compute the DFT for the averaged N samples.
 

Offline gf

  • Super Contributor
  • ***
  • Posts: 1114
  • Country: de
Let's assume that we are sampling N samples into a buffer. We have selected N in such way that ADC sampling frequency / LO frequency over N will become integer in all practical terms (for Fs=72Mhz/6/14=857142.8571428572_Hz and  LO=80kHz: Fs/LO=10.71428571428, we could select N = 93 * 10.71428571428 = 996.42857142857 = 996, for example). Now, we will compute the DFT for the given LO frequency of 80kHz (using quadrature correlator with pre-computed sin/cos tables with Hann-window). The RMS for this DFT = sqrt(sum(Q)^2 + sum(I)^2). No weird spectrum here.

Now, let's extend this a little, and that we are sampling K*N consecutive samples into a buffer. We will then divide this buffer into K blocks, each size of N samples. Then, we will compute DFT and RMS for each individual blocks of N samples. Finally, we will combine the results from all K blocks into one RMS value. There should not be no weird spectrum here, because we are just averaging K * individual DFT results. Right?

Vector averaging of the individual DFT results turns out to be mathematically equivalent to calculating the DFT over the whole K*N sized window, using a window function which is a K-fold repetition of the shorter N-tap window function. Just write-down the equation for sum that is calculated, for both cases, re-arrange them acordingly, and you'll see that both are the same. Consequence is that the equivalent filter frequency response is eventually that of the K-fold repeated N-tap window (i.e. the "weird" one).
Edit: I mean the DFT-term for the single frequency bin that is considered. The complete K*N sized DFT contains of course more frequency bins.

OTOH, if you calculate the RMS sum (power sum) of the individual DFT results (treating them as if the were independent and not necessarily coherent to each other), then the filter frequency response still corresponds to the (short) N-tap window function. The RMS sum does not contain phase information any more, and it reduces the variance of the noisefloor. It does not lower the noise floor. I.e. even if K approaches infinity, then N and the window function determine the lowest noise floor you can get.

Edit: Actually we'd like the get the filter response of a (non-repeated) window function, which spans the whole N*K window size, since its bandwidth is 1/K of the N-tap window's bandwidth then. But neither the vector nor the RMS averaging do that if the stored window function table has only N taps, and if the table is just indexed and not interpolated.

Edit: The ENBW of the "weird" frequency response is in fact as narrow as the ENBW of the desired frequency response. Unusual are the high side lobes. The main lobe is even narrower, conversely. If the signal spectrum happens to be "empty" (besides random noise) at the side lobe frequency (i.e. no spurs or harmonics), then this window function is still usable for the detector. It just requires a decent frequency planning to ensure that.

Btw, if you are willing to accept a gap at the center of the lobe, then I've a proposal for SA mode. Actually is is based on the above principle, with K=N=64, and doing RMS averaging of the individual sub-results. I've chosen the frequency plan to exclude the spurs and also the noise near DC. My expectation were a DR improvement of say 15+ dB, compared to simple RMS averaging over all samples. It is still rather wide-band (approx +/- 67kZh @-3dB), in order that the center gap remains small when compared to the total BW.
(Edit: The frequency plan is just an initial draft - there is still room for fine-tuning)

Edit: One more thing:
Quote
for Fs=72Mhz/6/14=857142.8571428572_Hz and  LO=80kHz: Fs/LO=10.71428571428, we could select N = 93 * 10.71428571428 = 996.42857142857 = 996, for example

In order that you can continue a vector sum coherently after accumulating N samples, by wrapping-around the sin/cos tables, an exact integral multiple of LO periods need to fint into the N samples. For fLO = 80kHz, this applies to N=75, or any integral multiple of 75, but not to 996.
« Last Edit: June 19, 2021, 02:07:22 pm by gf »
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
(Edit: The frequency plan is just an initial draft - there is still room for fine-tuning)

I have been thinking about this as well. Since the LO frequency needs to be multiple of 8000Hz, the current sampling frequency of 857142.857142 Hz is pretty unfriendly in this regard. The following options for the ADC sample rate are possible with different MCU core clock frequencies:

MCU core frequency 72MHz:
ADC sampling time 1.5cy: ADC sampling frequency 72MHz/6/14 = 857142.857142 Hz.
ADC sampling time 7.5cy: ADC sampling frequency 72MHz/6/20 = 600000Hz (75*8000Hz).

MCU core frequency 64MHz:
No good division ratio available which would provide advantage over 72MHz or 56MHz.

MCU core frequency 56MHz:
ADC sampling time 1.5cy: ADC sampling frequency 56MHz/4/14 = 1000000Hz (125*8000Hz).
ADC sampling time 7.5cy: ADC sampling frequency 56MHz/4/20 = 700000Hz (87.5*8000Hz).

Using the faster 72MHz core clock would provide faster computation times compared to 56MHz clock (kind of obvious). On the other hand, the 56MHz core clock would provide us faster ADC sampling rate and more samples in the given time period ie. better SNR.
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
Edit: One more thing:
Quote
for Fs=72Mhz/6/14=857142.8571428572_Hz and  LO=80kHz: Fs/LO=10.71428571428, we could select N = 93 * 10.71428571428 = 996.42857142857 = 996, for example

In order that you can continue a vector sum coherently after accumulating N samples, by wrapping-around the sin/cos tables, an exact integral multiple of LO periods need to fint into the N samples. For fLO = 80kHz, this applies to N=75, or any integral multiple of 75, but not to 996.

And this applies only using 72MHz and 7.5cy ADC sampling time (see my post above). If the system would use 72MHz and 1.5cy ADC sampling time, the 996 samples would not provide exact integer multiple of LO cycles. But it may be possible to nudge the ADC sample buffer pointer for each DFT computation so that the sampling will look practically coherent for each DFT computation.
 

Offline gf

  • Super Contributor
  • ***
  • Posts: 1114
  • Country: de
Edit: One more thing:
Quote
for Fs=72Mhz/6/14=857142.8571428572_Hz and  LO=80kHz: Fs/LO=10.71428571428, we could select N = 93 * 10.71428571428 = 996.42857142857 = 996, for example

In order that you can continue a vector sum coherently after accumulating N samples, by wrapping-around the sin/cos tables, an exact integral multiple of LO periods need to fint into the N samples. For fLO = 80kHz, this applies to N=75, or any integral multiple of 75, but not to 996.

And this applies only using 72MHz and 7.5cy ADC sampling time (see my post above). If the system would use 72MHz and 1.5cy ADC sampling time, the 996 samples would not provide exact integer multiple of LO cycles. But it may be possible to nudge the ADC sample buffer pointer for each DFT computation so that the sampling will look practically coherent for each DFT computation.

For NA usage, 12/14 MSa/s and 80kHz IF are fine. Seven 80kHz cyles == 75 samples then.

If a later SA goal were stiching the spectrum from multiple scans, then freqency steps which are in sync with DFT bins can be helpful.
Then 1MSa/s might be indeed be favorable, even if the CPU power is a bit lower then.

Edit: E.g. 4000 samples @1MSa/s would lead to "nice" 250Hz bins, and 8kHz were 32 bins then.
« Last Edit: June 19, 2021, 06:13:20 pm by gf »
 

Offline gf

  • Super Contributor
  • ***
  • Posts: 1114
  • Country: de
if you are willing to accept a gap at the center of the lobe, then I've a proposal for SA mode. Actually is is based on the above principle, with K=N=64, and doing RMS averaging of the individual sub-results. I've chosen the frequency plan to exclude the spurs and also the noise near DC. My expectation were a DR improvement of say 15+ dB, compared to simple RMS averaging over all samples. It is still rather wide-band (approx +/- 67kZh @-3dB), in order that the center gap remains small when compared to the total BW.

Attached is the expected frequency response of this RBW filter. As said, unfortunately with gap at the center :--. But there are too much spurs, noise and DC offset error in the center. Excluding it enables a higher dynamic range. IMO the gap needs to be filled eventually by stitching multiple measurements with frequency offset.

Edit: updated Link: https://godbolt.org/z/PG461P1jW
I expect a noise floor of roughly -67dBm (IF level), or a DR of almost 70dB if the max. IF level is 3dBm.

Edit: I played a bit with the gap width vs. DR trade-off.
At the cost of ~0.85dB the gap can be reduced by 40kHz, and reduction by 48 kHz costs ~1.5dB (see attached 2nd diagram).
(In the sin/cos tables this meens 36kHz or 32kHz instead of 56kHz)
« Last Edit: June 20, 2021, 07:45:16 am by gf »
 
The following users thanked this post: Kalvin

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
For NA usage, 12/14 MSa/s and 80kHz IF are fine. Seven 80kHz cyles == 75 samples then.

If a later SA goal were stiching the spectrum from multiple scans, then freqency steps which are in sync with DFT bins can be helpful.
Then 1MSa/s might be indeed be favorable, even if the CPU power is a bit lower then.

Edit: E.g. 4000 samples @1MSa/s would lead to "nice" 250Hz bins, and 8kHz were 32 bins then.

Good, let's keep the current sampling rate 72MHz/6/14.
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
if you are willing to accept a gap at the center of the lobe, then I've a proposal for SA mode. Actually is is based on the above principle, with K=N=64, and doing RMS averaging of the individual sub-results. I've chosen the frequency plan to exclude the spurs and also the noise near DC. My expectation were a DR improvement of say 15+ dB, compared to simple RMS averaging over all samples. It is still rather wide-band (approx +/- 67kZh @-3dB), in order that the center gap remains small when compared to the total BW.

Attached is the expected frequency response of this RBW filter. As said, unfortunately with gap at the center :--. But there are too much spurs, noise and DC offset error in the center. Excluding it enables a higher dynamic range. IMO the gap needs to be filled eventually by stitching multiple measurements with frequency offset.

Edit: updated Link: https://godbolt.org/z/PG461P1jW
I expect a noise floor of roughly -67dBm (IF level), or a DR of almost 70dB if the max. IF level is 3dBm.

Edit: I played a bit with the gap width vs. DR trade-off.
At the cost of ~0.85dB the gap can be reduced by 40kHz, and reduction by 48 kHz costs ~1.5dB (see attached 2nd diagram).
(In the sin/cos tables this meens 36kHz or 32kHz instead of 56kHz)

Thank you for investigating this. You are right that spurs and the frequency components near DC limit the dynamic range in SA mode. Filtering out the frequency components near DC means that there will be a gap in the spectrum near DC. I would see that allowing/tolerating the gap is a trade-off between better dynamic range and worse dynamic range in SA mode. It is quite possible to implement this as a feature that the user can select, because we can add new commands to serial port protocol quite easily. The compile-time configuration will then provide the default option* (either better dynamic range with a gap, or gapless spectrum with worse dynamic range).

Edit: I was playing with a simple M-tap comb-filter with null at DC, placed after IIR low-pass filter, which could be used as a DC-removal filter. The IIR low-pass filter would provide the adjustable RBW, and the comb-filter would take care of removing the low-frequency components near DC. Alternative solution would be to implement a IIR band-reject filter at DC with BW somewhere 2kHz ... 8kHz, but I have not investigated this. Lyons has some other ideas for narrow-band DC removal in his book "Understanding Digital Signal Processing", 3rd ed. section 13.40 "Linear-Phase DC Removal Filters" which shows some strategies for reducing the band-pass ripple of a typical comb-filter and making filter response narrower at DC.

Edit2: * In order to maintain compatibility with the existing PC software which do not support command protocol extansions, the compile-time configuration will be used as default.
« Last Edit: June 20, 2021, 09:33:06 am by Kalvin »
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
Vector averaging of the individual DFT results turns out to be mathematically equivalent to calculating the DFT over the whole K*N sized window, using a window function which is a K-fold repetition of the shorter N-tap window function. Just write-down the equation for sum that is calculated, for both cases, re-arrange them acordingly, and you'll see that both are the same. Consequence is that the equivalent filter frequency response is eventually that of the K-fold repeated N-tap window (i.e. the "weird" one).
Edit: I mean the DFT-term for the single frequency bin that is considered. The complete K*N sized DFT contains of course more frequency bins.

Ok, here is my attempt (see attachment). Let's assume that we are sampling our LO into a buffer of K*N samples, and the want to estimate the RMS of the periodic LO signal by splitting the effort into K individual blocks (each block has size N samples). Because we have chosen our ADC sample rate so that the N is always an integer multiple LO period (75 samples), we can combine individual one-bin DFT RMS values to get total signal RMS power. Due to a fact that we are actually performing coherent sampling, we can integrate the real and the imaginary parts separately over all blocks before computing the RMS of the signal, which will increase the SNR. But this will not introduce weird spectral response.
« Last Edit: June 20, 2021, 12:17:48 pm by Kalvin »
 

Offline gf

  • Super Contributor
  • ***
  • Posts: 1114
  • Country: de
Due to a fact that we are actually performing coherent sampling, we can integrate the real and the imaginary parts separately over all blocks before computing the RMS of the signal, which will increase the SNR. But this will not introduce weird spectral response.

In general this does introduce the weird response. Only for a rectangular window function it does not, since its K-fold repetition is still rectangular. See below.

Code: [Select]
// complex sine wave
lo = exp(-1i*2*pi*freq*[0:74]/fs);

// 4 chunks of samples
samp_1 = rand(1,75);
samp_2 = rand(1,75);
samp_3 = rand(1,75);
samp_4 = rand(1,75);

// regular coherent summation over all samples, note: window function spans all 75*4 samples!
w_all = hanning(75*4,"periodic")';
iq_all = sum([samp_1 samp_2 samp_3 samp_4] .* [lo lo lo lo] .* w_all);

// cohrerent summation in chunks
w = hanning(75,"periodic")';
iq_1 = sum(samp_1 .* lo .* w);
iq_2 = sum(samp_2 .* lo .* w);
iq_3 = sum(samp_3 .* lo .* w);
iq_4 = sum(samp_4 .* lo .* w);
iq_sum = iq_1 + iq_2 + iq_3 + iq_4;

// this is the same as
iq_sum = sum([samp_1 samp_2 samp_3 samp_4] .* [lo lo lo lo] .* [w w w w]);


But at the end iq_sum is not the same as iq_all !!!
Because [w w w w] != w_all.
While fft(w) and fft(w_all) look as usual, fft([w w w w]) is the "weird" frequency response, with the extra side lobe.

OTOH, the following would lead to iq_sum == iq_all:
Code: [Select]
w_1 = w_all(1:75);
w_2 = w_all(76:150);
w_3 = w_all(151:225);
w_4 = w_all(226:300);

iq_1 = sum(samp_1 .* lo .* w_1);
iq_2 = sum(samp_2 .* lo .* w_2);
iq_3 = sum(samp_3 .* lo .* w_3);
iq_4 = sum(samp_4 .* lo .* w_4);
iq_sum = iq_1 + iq_2 + iq_3 + iq_4;

But this requires to store the whole N*K sized window function table, in order that w_1 ... w_k can be extracted.
« Last Edit: June 20, 2021, 02:28:43 pm by gf »
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
I think we have had some misunderstanding, or I have not been able to explain my intention properly. However, you have captured my intention here in compact format:

Code: [Select]
// this is the same as
iq_sum = sum([samp_1 samp_2 samp_3 samp_4] .* [lo lo lo lo] .* [w w w w]);

If the system has a maximum sample buffer size of K*N samples, then the user may select to capture (1, 2, 3 or ... K) * N samples, depending how fast sweep step rate and dynamic range the user wants to have. Thus the system is required to store only on set of windowed sin and windowed cos tables of N samples for all different sweep rates.

I was originally thinking that the system would provide four different sweep rates, for example 0.5K, 1K, 2K and 4K samples, and the system would need to have four set of windowed sin/cos tables: 0.5K, 1K, 2K and 4K in size. Since this is not very flexible and would waste Flash memory because of four set of tables, I thought that the K*N sample approach would be more flexible and more economic.

Edit: This is something that I do not understand again: While fft(w) and fft(w_all) look as usual, fft([w w w w]) is the "weird" frequency response, with the extra side lobe. Why this fft([w w w w]? At least my intention is more like fft(w)+fft(w)+fft(w)+fft(w), like averaging four different FFTs/DFTs. Since we are capturing in coherent way, we can integrate over all four FFTs/DFTs in order to get better dynamic range, compared to four individual FFTs/DFTs.
« Last Edit: June 20, 2021, 05:22:31 pm by Kalvin »
 

Offline gf

  • Super Contributor
  • ***
  • Posts: 1114
  • Country: de
I was originally thinking that the system would provide four different sweep rates, for example 0.5K, 1K, 2K and 4K samples, and the system would need to have four set of windowed sin/cos tables: 0.5K, 1K, 2K and 4K in size. Since this is not very flexible and would waste Flash memory because of four set of tables, I thought that the K*N sample approach would be more flexible and more economic.

Yes, I understand the reason, but

the result having the intended frequency response were rather this one:
Code: [Select]
N=75; K=4;
w_all = hanning(N*K,"periodic")';
iq_all = sum([samp_1 samp_2 samp_3 samp_4] .* [lo lo lo lo] .* w_all);

while the following calculatioin leads to the "weird" frequency response with the extra side lobes:
Code: [Select]
iq_sum = sum([samp_1 samp_2 samp_3 samp_4] .* [lo lo lo lo] .* [w w w w]);

I don't say that the "weird" freqency response won't work. It still can work, too.
But the latter filter is not equivalent to the former one.


Quote
This is something that I do not understand again: While fft(w) and fft(w_all) look as usual, fft([w w w w]) is the "weird" frequency response, with the extra side lobe. Why this fft([w w w w]? At least my intention is more like fft(w)+fft(w)+fft(w)+fft(w), like averaging four different FFTs/DFTs. Since we are capturing in coherent way, we can integrate over all four FFTs/DFTs in order to get better dynamic range, compared to four individual FFTs/DFTs

General form of single DFT term (with window function) is basically
Code: [Select]
sum(SAMPLES .* LO .* WINDOW)
where SAMPLES, LO, WINDOW are row vectors of equal size (-> number of samples).
LO is the complex sine wave. The frequency response of the corresponding filter is fft(WINDOW) 1).

When you calculate
Code: [Select]
iq_1 = sum(samp_1 .* lo .* w);
iq_2 = sum(samp_2 .* lo .* w);
iq_3 = sum(samp_3 .* lo .* w);
iq_4 = sum(samp_4 .* lo .* w);
iq_sum = iq_1 + iq_2 + iq_3 + iq_4;
which is the same as
Code: [Select]
iq_sum = sum([samp_1 samp_2 samp_3 samp_4] .* [lo lo lo lo] .* [w w w w]);
then you effectively calculate a single DFT term for SAMPLES=[samp_1 samp_2 samp_3 samp_4], LO=[lo lo lo lo] and WINDOW=[w w w w]
I.e. the effecive window function for this single-term DFT calculation is [w w w w] then.
Therefore the frequency response of the corresponding filter is fft([w w w w]).

When you do the following "brute force" calculation where you sweep frequency and calculate the vector-averaged detector response for each frequency point, then you get the same frequency response as from fft([w w w w]). So I think it can't be so wrong.
Code: [Select]
fs = 12000000/14;
freqs = 30000:100:130000;
lo = exp(-1i*2*pi*80000*[0:74]/fs);
w = hanning(75,"periodic")';
w /= sum(w);

H = freqs .* 0;
for i = 1:length(freqs)
  f = freqs(i);
  samples = sin(2*pi*f*[0:299]/fs);
  for j=0:75:299
    H(i) += sum(samples(j+1:j+75) .* lo .* w) / 2;
  end
end

clf
plot(freqs,20*log10(abs(H)));
ylim([-100 0])
xlim([freqs(1) freqs(end)])
grid on


1) Plot frequency responses. Zero-fill to align different lengths, and to interpolate for a higher freqency resolution.
Code: [Select]

fs = 12000000/14;

% approx. 100Hz resolution for the plot
NFREQ=floor(fs / 100 + 0.5);            % number of frequency points
f = [0:NFREQ-1] / NFREQ * fs - fs / 2;  % frequency scale

% the window functions we want to plot
w = hanning(75, "periodic")';
w_all = hanning(75*4, "periodic")';
w4 = [w w w w];

% normalize to sum==1, in order that we get 0dB for DC
w /= sum(w);
w_all /= sum(w_all);
w4 /= sum(w4);

% zero-fill in order to interpolate and align different lengths
w = [w zeros(1,NFREQ-length(w))];
w_all = [w_all zeros(1,NFREQ-length(w_all))];
w4 = [w4 zeros(1,NFREQ-length(w4))];

clf; hold on;
plot(f,20*log10(abs(fftshift(fft(w)))),";w;")
plot(f,20*log10(abs(fftshift(fft(w_all)))),";w-all;")
plot(f,20*log10(abs(fftshift(fft(w4)))),";w4;")
xlim([-50000 50000])
ylim([-100 0])
grid on

« Last Edit: June 21, 2021, 07:11:04 am by gf »
 

Offline gf

  • Super Contributor
  • ***
  • Posts: 1114
  • Country: de
Kalvin, here is yet another idea to calculate the 75*K point single-frequency DFT term with Hanning window, for any K > 1, using only 900 bytes of table data.
Note that It is also possible to apply the window function in the frequency domain, after performing the DFT.
fft(a.*b) is equivalent to circular_convolution(fft(a), fft(b)), consequently fft(samples.*window) == circular_convolution(fft(samples), fft(window)).
The DFT of a periodic Hanning window fft(hanning(M,"periodic")')/M is always [ 0.5 -0.25 0 ... 0 -0.25 ], for any M >= 3.
This means the Hanning window can easily be applied in the frequency domain, by (circular) convolution with [ -0.25 0.5 -0.25 ]. Fortunately, these are only 3 non-zero taps :)
In order that this convolution can be calculated at the desired frequency, the DFT term needs to be calculated not only for the desired frequency, but also for the lower and upper neighbor bins of the desired frequency.
This is computationally more expensive (a couple of additional instructions per sample), but it saves a significant amount of memory.
The tables need to be calculated individually for each K, so they cannot reside in ROM, but must be in RAM. OTOH, only 75*6*2 = 900 bytes are required.


Edit:

I have to retract my statement partially. A window function can indeed be applied in the frequency domain as well. But unfortunately it does not work out with the sizes as I was thinking. Sorry for my calculation mistake. While the basic priciple is still valid, a phase correction seems unavoidable for each chunk, which comes at the cost of two additional complex multiplications per chunk. My aim is still to keep the algorithm so cheap that it can be calculated in real-time.
« Last Edit: July 04, 2021, 11:19:48 am by gf »
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
gf, thank you for your ideas!  :-+ Sorry that I have not been active for a week or so, but I have now time to proceed with this little project again.

I just added the documentation for the hardware modifications at Github:
https://github.com/kalvin2021/ltdz-dsp

The documentation for building the actual firmware is still in progress.



 

Offline gf

  • Super Contributor
  • ***
  • Posts: 1114
  • Country: de
Sorry for the mistake in my previous statement which I did retract.
I'm still waiting for my order to arrive. Likely I won't have much free time during summer either to play with it when it arrives.
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
gf, You asked whether the LTDZ board could be used for measuring narrow band crystal filters. I read the ADF4351 datasheet and wrote the following script for testing the register configurations. It seems that it is possible in theory to generate frequency steps in 10Hz resolution (or even better), but I have not tested this in practice. Hopefully I didn't miss something.

Code: [Select]
# Reference oscillator frequency
global REF_Hz = 25000000;

# Reference doubler: 0|1
global D = 0;

# Refecence divide by 2: 0|1
global T = 1;

# Preset divide ratio in range [1, 1023]
global R = 100;
 
# Modulus in range [2, 4095]
global MOD = 4095; 

# PFD frequency
global pfd_Hz = REF_Hz * ((1 + D)/(R * (1 + T)));

# Compute ADF4351 INT and FRAC register values for the given frequency.
#
# f       Desired frequency 35MHz <= f < 4400MHz.
# INT     INT-register value.
# FRAC    FRAC-register value.
# rf_Hz   Actual output frequency.
# err_Hz  Frequency error.
#
function [INT, FRAC, rf_Hz, err_Hz] = adf4351_reg(f)
  global REF_Hz;
  global MOD;
  global D;
  global T;
  global R;
  global pfd_Hz;
 
  if f <= 68750000
    RFDIV = 64;
  elseif f <= 137500000
    RFDIV = 32;
  elseif f <= 275000000
    RFDIV = 16;
  elseif f <= 550000000
    RFDIV = 8;
  elseif f <= 1100000000
    RFDIV = 4;
  elseif f <= 2200000000
    RFDIV = 2;
  else
    RFDIV = 1;
  endif

  VCO = f * RFDIV;
  assert(2200000000 <= VCO && VCO <= 4400000000);
 
  INT  = fix(VCO / pfd_Hz);
  assert(75<= INT && INT <= 65535);

  FRAC = round(((VCO / pfd_Hz) - INT) * MOD);
  assert(0 <= FRAC && FRAC <= (MOD-1));

  rf_Hz = round(REF_Hz * ((1 + D)/(R * (1 + T))) * (INT + FRAC / MOD) / RFDIV);

  err_Hz = f - rf_Hz;
endfunction

For example:
Code: [Select]
>> [INT, FRAC, f, err] = adf4351_reg(35000000)
INT =  17920
FRAC = 0
f =  35000000
err = 0
>> [INT, FRAC, f, err] = adf4351_reg(35000001)
INT =  17920
FRAC =  2
f =  35000001
err = 0
>> [INT, FRAC, f, err] = adf4351_reg(35000002)
INT =  17920
FRAC =  4
f =  35000002
err = 0
« Last Edit: July 04, 2021, 01:46:18 pm by Kalvin »
 

Offline gf

  • Super Contributor
  • ***
  • Posts: 1114
  • Country: de
R=100 is quite much. I wonder what's the phase noise then. My understanding is: Lower fPFD -> larger INT -> more phase noise. To mitigate phase noise then, I also wonder whether the loop filter bandwidth would need to be reduced (significantly), at the expense of a larger lock time.
See also https://www.analog.com/media/en/training-seminars/tutorials/MT-086.pdf
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
R=100 is quite much. I wonder what's the phase noise then. My understanding is: Lower fPFD -> larger INT -> more phase noise. To mitigate phase noise then, I also wonder whether the loop filter bandwidth would need to be reduced (significantly), at the expense of a larger lock time.
See also https://www.analog.com/media/en/training-seminars/tutorials/MT-086.pdf

Good points! I did not consider the performance at this point, so I overlooked that issue. It is possible to experiment with some other R values in range [10, 150] as well. The point here is that it is quite possible to perform sweeps with smaller steps. For limited frequency range(s), the R can be hand-selected so that the performance will be optimized in that regard.
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
Added a simple Python/Tkinter-utility for reading and plotting LTDZ ADC sample buffer data. This is a development tool for the project development purposes, not intended to be used as a general purpose spectrum analyzer or network analyzer tool. The sweep data can be exported to MATLAB and GNU/Octave for later analysis.

Project page:
https://github.com/kalvin2021/ltdz-dsp

Edit: Changed Matplotlib style/theme.
« Last Edit: July 11, 2021, 01:39:53 pm by Kalvin »
 

Offline KalvinTopic starter

  • Super Contributor
  • ***
  • Posts: 2145
  • Country: fi
  • Embedded SW/HW.
Added the firmware build instructions in the project page at Github:
https://github.com/kalvin2021/ltdz-dsp
 
The following users thanked this post: gf, mahdiks


Share me

Digg  Facebook  SlashDot  Delicious  Technorati  Twitter  Google  Yahoo
Smf