Technique of sampling a high-frequency signal (Mhz) with a lower sampling rate (khz)?

Question

I currently have a .wav signal that was recorded at a 48 kHz sample rate, with a central frequency of 5.260 MHz and bandwidth of 4 kHz. I'm trying to get some information from this signal, specifically the location of the frequency peaks using matlab, which I would expect to be around 5.260 MHz. The code I'm using to generate the power spectral density is as follows:

path = '5.260MHz.wav';
f0 = 5.260e6; % frequency 5.260 MHz

% sample properties
[x, Fs] = wavread(path); % Fs = 48 kHz
alias = 48e3*floor(f0/48e3);      % starting frequency of target alias

%estimate spectrum
[psd, f] = pwelch(x, 512, [], [], Fs);

% If the frequency lies in one of the mirror frequency bands,
% we have to rotate about half the sample rate.
if (mod(f0, 48e3) > 24e3)  
    f = Fs - f;
end

plot(f + alias, psd)

This script produces the following plot

The bandwidth does seem to be 4 kHz as expected, but the frequency peaks appear to be in the wrong place. I would expect them to occur between 5.256 and 5.264 MHz, or maybe 5.258 and 5.262 MHz, not between 5.276 and 5.280 MHz.

However, when I generated my own files with the same central frequency and bandwidth, I got the following when I ran the above script:

% signal properties
f0 = 5.260e6;           % frequency 5.260 MHz
sf = 4e3;               % bandwidth 4 kHz

% sample properties
fs = 24e6;              % sampling frequency 24 MHz
N = 10e6;               % number of samples

% generate random frequency modulated sinusoidal signal
fi = smooth(randn(N, 1), 11);
fi = fi / std(fi) * sf + f0;
x = cos(2 * pi * cumsum(fi) / fs);

% downsample to 48 kHz (factor 500)
x = x(1 :500: end);
fs = fs / 500;

wavwrite(x, fs, 'test_5.260MHz.wav');

This plot is exactly as expected, using the same script as the first plot. I'm wondering if I should be treating the first input different in some way to have it graphed correctly.

So far I've tried testing some other generated cosine signals using the same code as above, and they all graph correctly. However, every data signal I've tried has been off to some degree (I can post different sample plots if that would be useful). It could be an error with the data collection instrument, but it's far more likely to be a bug in my code.

Answer 1

You cannot capture a 5 MHz signal by sampling at 48 kHz. That doesn't mean there won't be anything in your sample, it means it will represent the signal incorrectly.

What happens if you use a too small sampling frequency is called aliasing: the original frequency shows up as another frequency altogether. The translation from the real to the sampled frequency is like a process of wrapping or folding the infinite frequency axis into your limited frequency interval $[0 ; \frac12 f_s]$, or $[-\frac12 f_s ; \frac12 f_s]$ for complex-valued signals. $\frac12 f_s$ is also called Nyquist frequency.

In your example, not just the main frequency but the total bandwidth of the signal of 800 kHz is much larger than the sampling frequency. Because of this, I wouldn't expect any clearly defined peak to be present. Since you seem to find peaks, my guess is that that is due to your very primitive approach to spectral estimation. If you use a better method, e.g. Welch's modified periodogram, which is implemented in Matlab as pwelch, I'd expect that you see a very flat (white) spectrum.

---

Addition based on nimrodm's comment: My initial statement has to be relativized if the bandwidth of the signal is smaller than $\frac12 f_s$. In that case, if there is additional knowledge about the true signal frequency, this knowledge can be used to resolve the ambiguity introduced by aliasing, and perfect reconstruction of the true signal is possible. See undersampling. However, it is not possible to determine the true frequency from the undersampled signal and nothing else.

---

A simulated example: Matlab code to generate a frequency-modulated signal with peak at 5 MHz and bandwidth 400 kHz, initially sampled at 24 MHz.

% signal properties
f0 = 5e6;               % frequency 5 MHz
sf = 0.4e6;             % bandwidth 400 kHz

% sample properties
fs = 24e6;              % sampling frequency 24 MHz
N = 10e6;               % number of samples

% generate random frequency modulated sinusoidal signal
fi = smooth(randn(N, 1), 11);
fi = fi / std(fi) * sf + f0;
x = cos(2 * pi * cumsum(fi) / fs);

Estimating the spectrum using pwelch

% estimate spectrum
[psd, f] = pwelch(x, 512, [], [], fs);
figure
plot(f, psd)
xlabel('f [Hz]')
ylabel('normalized power spectral density [1/Hz]')

gives this result:

If the signal is downsampled to 48 kHz

% downsample to 48 kHz (factor 500)
fs = fs / 500;
x = x(1 :500: end);

the estimated spectrum looks like this:

This is, as expected, basically a white (flat) spectrum, and the many small "peaks" are random fluctuations due to the fact that this spectrum was estimated based on only a 500th of the original 10 million samples. The estimate could be improved by repeating the procedure and averaging over the different results.

If the bandwidth is changed to 4 kHz

sf = 4e3;               % bandwidth 4 kHz

the results of spectal estimation on original and downsampled signal look like this:

The narrower peak fits nicely into the range 0 to 24 kHz of the downsampled spectral estimate, it just necessarily appears at the wrong place, in this case 8 kHz.

8 kHz = 5 MHz – 104 · 48 KHz is the smallest alias of 5 MHz with respect to a sampling rate of 48 kHz, the next largest is 48 kHz – 8 kHz = 40 kHz (mirroring around the Nyquist frequency of 24 kHz). All the larger aliases derive from these two by adding integer multiples of 48 kHz: 56, 88, 104, 136 kHz up to 4952, 4984, 5000, 5032 kHz and so on.

Since the spectral estimation cannot distinguish these aliases, the frequency axis of the spectrum plot becomes ambiguous: instead of representing frequency components of the signal from 0 to 24 kHz, it could as well represent frequencies from 48 to 24 kHz (backwards because of the mirroring!),

or from 48 to 72 kHz, and so on.

If we know the frequency content of the signal lies around 5 MHz, we can select the corresponding aliased frequency interval, 4992 to 5016 kHz,

and recover the correct location of the peak at 5 GHz.

This trick does not work, however, if the bandwidth of the signal is larger than $\frac12 f_s$, because the true frequency content of the signal has been aliased across several of these frequency intervals.

Answer 2

You can capture a higher frequency signal by a using a sample rate much lower than twice the signal/spectral content frequency. It's called undersampling. To be useful, one has to know approximately how many "folds" the signal frequency is above the sample rate, there has to be no baseband (or other folded spectrum) signal to alias against, thus the bandwidth of the signal has to be narrower than Fs/2, and the spectrum of the signal can't cross any multiple of Fs/2. The jitter of the sampling clock also has to be quite low (as low as the max absolute time jitter required for a sampler that is clocked at well above the signal frequency.)

This is no different from baseband sampling, as one needs a-priori knowledge (that the signal was, in fact, correctly low-pass-filtered baseband below Fs/2) to reconstruct the original signal correctly,