Oversampling and decimation leads to slightly different result in frequency domain

Question

I am trying to explore the effect of oversampling and decimation with the following Matlab code

% Step 1: generate two sequences of data (here, I consider them to be frequency-domain data) with only the first element being the same
a = randn(36,1) + 1j*randn(36,1); 
b = randn(36,1) + 1j*randn(36,1); 
b(1) = a(1) 
% Step 2: upsample and convert to the time-domain
a1 = ifft(a,36*6);
b1 = ifft(b,36*6);
% Step 3: convert the time-domain to frequency domain. As expected, we have a2(1) = b2(1) = a(1) = b(1)
a2 = fft(a1);
b2 = fft(b1);
% Step 4: decimate  (filtering and downsampling) the time-domain samples by a factor and then convert to frequency domain again
a3 = decimate(a1,6); b3 = decimate(b1,6);
a4 = fft(a3); b4 = fft(b3)

Now, we can see that $a4(1) \neq b4(1)$, this is weird as both frequency domain sequences have undergone the filters with the exactly the same frequency responses, so the effect on the data at the same frequency position should be the same. However, this inequality shows that this is not the case.

Also, we can see that $a4(1)$ or $b4(1)$ is not equal to $a(1) = b(1)$ any more.

I have difficulty with understanding these two things. Hope someone could help me with understanding this. Thanks.

Answer 1

Oversampling and decimation leads to slightly different result in frequency domain

Which is expected. As long as decimation involves low-pass filtering, the low-pass filter will change the signals to some extent. The only way around this would be an ideal lowpass filter, which is infinitely non-causal so it cannot be implemented in the real world.

Answer 2

There are two aspects that need to be corrected: (1) upsampling in the frequency domain, which needs to be done in consideration of the fft frequencies and scaling; and (2) you must not low pass filter prior to decimation because doing so is going to affect some of your frequency content (note that you can get away with this because upsampling in the frequency domain does not introduce any frequency content above the original Nyquist frequency, thus there is no risk of aliasing when you downsample back to the original sample rate).

In general, be aware that the fft result is proportional to the signal length, thus your comment As expected, we have a2(1) = b2(1) = a(1) = b(1) is actually not at all expected if the upsampling was done correctly, so was one of the indications that something was incorrect.

Here is corrected Matlab code that solves your problems thereby obtaining the expected result that upsampling in the frequency domain followed by downsampling in the time domain recovers the original signal:

% Step 1: generate two sequences of data (here, I consider them to be frequency-domain data) with only the first element being the same
n = 36;
a = randn(n,1) + 1j*randn(n,1);
b = randn(n,1) + 1j*randn(n,1);
b(1) = a(1);

% Revised Step 2: upsample and convert to the time-domain
up_sample_factor = 6;   % must be an integer > 1
nyquist_index = n/2+1;
if rem(n,2)
    % length is odd, so it does not include the nyquist frequency
    tmp_a = [a(1:floor(nyquist_index)); zeros(n*(up_sample_factor-1),1); a(ceil(nyquist_index):end)]*up_sample_factor;
    tmp_b = [b(1:floor(nyquist_index)); zeros(n*(up_sample_factor-1),1); b(ceil(nyquist_index):end)]*up_sample_factor;
else
    % length is even, so it includes the nyquist frequency
    tmp_a = [a(1:nyquist_index-1); a(nyquist_index)/2; zeros(n*(up_sample_factor-1)-1,1); a(nyquist_index)/2; a(nyquist_index+1:end)]*up_sample_factor;
    tmp_b = [b(1:nyquist_index-1); b(nyquist_index)/2; zeros(n*(up_sample_factor-1)-1,1); b(nyquist_index)/2; b(nyquist_index+1:end)]*up_sample_factor;
end
a1 = ifft(tmp_a);
b1 = ifft(tmp_b);

% Step 3: convert the time-domain to frequency domain. As expected, we have a2(1) = b2(1) = up_sample_factor*a(1) = up_sample_factor*b(1)
a2 = fft(a1);
b2 = fft(b1);

% Revised Step 4: decimate (without low pass filtering) the time-domain samples by a factor and then convert to frequency domain again
a3 = downsample(a1,up_sample_factor);
b3 = downsample(b1,up_sample_factor);
a4 = fft(a3); b4 = fft(b3);

% final check
ifft_a = ifft(a);
ifft_b = ifft(b);
[ifft_a, a3]    % these are now correctly identical
[ifft_b, b3]    % these are now correctly identical
[a(1); b(1); a2(1)/up_sample_factor; b2(1)/up_sample_factor; a4(1); b4(1)]  % these are now correctly identical