Additional peak after resampling

Question

I prepared a Matlab script, which re-samples spectrum from 24kHz to 20kHz. When I used up-sampling, a fir filter, and then up-sampling again at the end, I get an additional peak at 4kHz which is not expected. When I checked via the resample function from Matlab paltette, there is no additional peak. What did I do wrong? What do I have to add to my script to avoid this kind of situation. Thank you in advance

clc;
close all;
clear variables;
Fs1 = 24000;
TimeStep1 = 1/Fs1;
Lenght = 10000;
time1 = (0:Lenght-1)*TimeStep1;
InputSignal = sin(2*pi*8000*time1);
NFFT1 = 2^nextpow2(Lenght);
Y1 = fft(InputSignal,NFFT1)/Lenght;
f1 = Fs1/2*linspace(0,1,NFFT1/2+1);
figure
plot(f1,2*abs(Y1(1:NFFT1/2+1))) 
title('Original Spectrum')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
Fs2 = 20000;
T2 = 1/Fs2;
nGcd = gcd(Fs1, Fs2);
upY = upsample(InputSignal,Fs2/nGcd);
firCoeff = fir1(255, 0.9, 'low');
filteredY = filter(firCoeff,1,upY);
dataOut = resample(InputSignal,Fs2,ceil(Fs1));
L2 = length(dataOut);
t2 = (0:L2-1)*T2;
N=2;
lenght = (L2/N);
NFFT = 2^nextpow2(lenght);
Y = fft(dataOut,NFFT)/lenght;
f = Fs2/2*linspace(0,1,NFFT/2+1);
figure
plot(f,2*abs(Y(1:NFFT/2+1))) 
title('Resampled Spectrum')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
figure
dataOut2 = downsample(filteredY,Fs1/nGcd);
Y2 = fft(dataOut2,NFFT)/lenght;
plot(f,2*abs(Y2(1:NFFT/2+1)))
title('UpFIRDwn Spectrum')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

Answer 1

You are upsampling by a factor 5, which creates mirrors of the original spectrum (see the spectrum of $Yup$). With the filter you intent to remove these mirrors such that you can safely downsample the signal without aliasing.

To prevent aliasing, the upper limit for your cutoff frequency is half the lowest sampling rate, which is $0.5*Fs2 = 10 kHz$. However, the finite steepness of your filter creates a so-called transition band, which leads to some aliasing for higher frequencies for this maximum cutoff frequency. Therefore, a cutoff frequency a little below the absolute maximum is desired, for example 9 kHz. Choosing the transition band depends on your filter length and singal content.

For a cutoff frequency of 9 kHz at 5 * 24 / 2 khz, the normalized cut-off frequency is 0.15. You can use $Fs2/nGcd * Fs2$ to calculate this normalized cutoff frequency automatically.

Some guessing from my side: In your code you use a normalized cutoff frequency of 0.9, which might be related to the 9 kHz and the output sampling rate of 20 kHz. However, the filter is applied to the $Yup$ signal, which has a sampling rate of 5 * 24 kHz.

In short: changing the normalized cutoff frequency of your filter from 0.9 to 0.15 should solve your problem.

EDIT: Follow up on using fir2

With the fir2(N, F, A) method of Matlab you can design an N points FIR filter by prescribing the filter magnitude at certain frequencies. In your case, you want to design a low pass filter, which means that the frequencies and magnitudes should be something like the following:

A = [1 1 0 0];
F = [0 wpass wstop 1];

where wpass and wstop are the relative frequencies up to where you want the filter to let the signal pass or stop, respectively.

In your case you can set these for example as follows:

wstop = Fs2/(Fs2/nGcd*Fs1);
wpass = 0.9 * wstop;

Selecting the filter length N is a tradeoff between complexity and the ripple of the filter with respect to an ideal low pass filter.

Note: In your question you have m=[0 0 1 1], which ticks on a high pass filter.