Separating two sound file frequencies using bandpass filter matlab

can someone please help me? I have a project to implement a superheterodyne transmitter and receiver without using matlab functions that modulate. The first part which I can't seem to get right is the RF filter, which the my doctor said should just be a bandpass filter centered around Wc. I am supposed to get two sound signals, add them together which is my TdmSig then pass them through the bandpass filter centered around the frequency of the first carrier signal (which modulates signal 1) and should remove signal 2.

The center of my bandpass filter is 100khz/2pi = Fc = 15900

My Fs1 = Fs2.. and Fs1 = 44100 before interp command and after Fs1 = 132300

I chose a bandwidth of 4000. I really tried for 3 days and I dont know how to get the bandwidth of the signal.

Also before all that I am worried maybe I might be modulating my signals incorrectly I dont understand very well

here is file he sent

The AM modulator: the modulation type is double side band suppressed carrier (DSB-SC) for all the signals, that is the modulation is performed by a simple multiplication with a carrier sinusoid. The first signal is modulated with 100KHz carrier. Each following signal is modulated with a carrier of frequency wc = 100khz + n*df where df= 60 KHz and the index n is the signal index (n = 0 for the first signal that is modulated at 100 KHz). The modulated signals are used to form a TDM signal

The RF stage: As discussed in the lecture, this is the stage that performs interference-image rejection. For simplicity, this stage will be implemented as a band-pass Filter (BPF) only, centered at the carrier frequency (that is tunable at the desired station). You can refer to MathWorks tutorial for the design of BPF in MATLAB:

my code

close all
clear
clc
[stereo1 Fs1] = wavread('Short_FM9090'); % This will save the samples of the 2 channels of the signal in the array stereo
[stereo2 Fs2] = wavread('Short_RussianVoice');

% Part i a(Converting to mono)
mono1(:,1)= stereo1(:,1)+stereo1(:,2); %This will add the 2 channels frequencies and convert to mono
mono2(:,1)= stereo2(:,1)+stereo2(:,2);

%take small part for fast computation
mono1 = mono1(1:200000);
mono2 = mono2(1:200000);

l1 = length(mono1); %This will return the size of the array mono1
l2 = length(mono2);

%Part ii b Amplitude Modulation

Mod1 = interp(mono1,3);  % increase samples for Nyquist to be satisfied
Mod2 = interp(mono2,3);
Fs1 = Fs1*3;
l1 = l1*3;

df = 60*10^3;
wc = 100*10^3;

n = 0:0.05:l1/20;
n=n(1,1:end-1);

carrier1 =  cos(wc*n/Fs1)';
Mod1 = Mod1.*carrier1;

%plot carrier

% figure
% plot(carrier1)
% title('carrier 1')

Mod1_fft = fft(Mod1);

zdouble = Mod1_fft;
zdouble = abs(zdouble(1:length(zdouble)/2+1));
frqdouble = (0:length(zdouble)-1)*Fs1/length(zdouble)/2;

subplot(2,1,1); plot(frqdouble,zdouble);
title('Spectrum of double-sideband Mod 1 signal');

wc2 = wc+df;

carrier2 = cos(wc2*n/Fs1)';
Mod2 =  Mod2.*carrier2;

Mod2_fft = fft(Mod2);
zdouble = Mod2_fft;
zdouble = abs(zdouble(1:length(zdouble)/2+1));
frqdouble = (0:length(zdouble)-1)*Fs1/length(zdouble)/2;

subplot(2,1,2); plot(frqdouble,zdouble);
title('Spectrum of double-sideband Mod 2 signal');

figure
subplot(2,1,1); plot(Mod1);
title(' Mod 1 signal');
subplot(2,1,2); plot(Mod2);
title('Mod 1 signal');

TdmSig =Mod1+Mod2;

figure
plot(TdmSig)
title( 'TDM Signal in Time domain before filter')
ylabel('Amplitude')
xlabel('Time')

Tdm_fft = fft(TdmSig);
Tdm_fft = fftshift(Tdm_fft);
k = -l1/2:(l1/2)-1;
figure
plot(k,abs(Tdm_fft))
title('Tdm in frequency domain before filter')

%Part iv (RF Filter centered at Wc and Tunable)

FsFilter =Fs1; % Sampling frequency
A_stop1 = 60; % Attenuation in the first stopband
F_stop1 = 11000; % Edge of the stopband
F_pass1 = 12900;% Edge of the passband
F_pass2 = 16900;% Closing edge of the passband
F_stop2 = 20000;% Edge of the second stopband
A_stop2 = 60; % Attenuation in the second stopband
A_pass = 1; % Amount of ripple allowed in the passband

BandPassSecObj = fdesign.bandpass(F_stop1,F_pass1,F_pass2,F_stop2,A_stop1,A_pass,A_stop2,FsFilter);
hd = design(BandPassSecObj,'butter');
demod= filter(hd,TdmSig);

soundsc(TdmSig) % should at least remove unwanted sound file

figure
plot(demod)
title('Signal after Bandpass filter for TDM signal')

Tdm_fft = fft(demod);
Tdm_fft = fftshift(Tdm_fft);
k = -l1/2:(l1/2)-1;
figure
plot(k,abs(Tdm_fft))
title('Tdm in frequency domain after filter')


I did not follow/check all of your code, but I have the impression that there is some confusion about frequency, angular frequency, the place of the factor $2\pi$ and the arguments of the trig functions.
For a frequency $f$ in cycles per second or Hertz, one possible signal exhibiting this frequency is $\cos(2\pi f\cdot t)$ so that $ω=2\pi f$ is the angular frequency in radians per second.
Now you want to shift the frequency of a base band signal $x$ to a frequency band with carrier frequency $f_c$ and thus angular frequency $ω_c=2\pi f_c$, which is accomplished by computing the product $y(t)=x(t)\cos(2\pi f_c\cdot t)$. Which gives in fact two copies of the basis signal, one around $f_c$ and one around the mirror frequency $-f_c$.
Since $x$ is sampled with sampling frequency $f_s$, and thus period $dt=\tfrac1{f_s}$, $x[n]=x(n/f_s)$, the samples of $y$ will be $y[n]=x[n]\cdot \cos(2\pi f_c/f_s\cdot n)$. Or with $ω_c/f_s$ instead of $2\pi f_c/f_s$, $y[n]=x[n]\cdot \cos(ω_c/f_s\cdot n)$
It is not comprehensible why the number $1/20$ as step size in $n$ occurs in this context, $n$ is here an integer, i.e., with step size $1$.