# Oversampling and decimation leads to slightly different result in frequency domain

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.

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.

• But, I expect the nonideal low-pass filtering is the same for both data, which means that the effect on a4(1) and b4(1) should be the same since they are at the same frequency position.
– Vic
Commented Jun 18 at 2:32
• Mostly truncation effects. Your signal is very short, so the filter tail that will get cut off will result in a significant error which is signal dependent. Another reason is that the signal is too short for the filter to do it's job, so you will end up with significant aliasing that's also signal dependent. Both effects will affect the first frequency Commented Jun 18 at 21:21

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

• I’m a little worried a beginner or student could read this as “do not low pass filter before decimation”, which is of course incorrect in most cases. I think the OP was wondering why he was seeing the discrepancy when seemingly applying the same routine to two different signals. I’m not sure “fixing” that discrepancy is of any value, rather the explanation is.
– Jdip
Commented Jun 20 at 6:50
• To address that concern, I've added context to that statement and also stated the conditions under which it is ok. Commented Jun 20 at 7:26