Hi Abdullah,
A couple of issues to address in the anlaysis.
I think you started by trying to define an ideal filter in the frequency domain based on its Continuous Time Fourier Transform (CTFT), like so.
syms w t real
H(w) = rectangularPulse(2*sym(pi)*720,2*sym(pi)*740,w);
Its impulse response is
h(t) = rewrite(ifourier(H(w),w,t),'exp');
[num,den] = numden(h(t));
h(t) = num/den
The code in the question defined h(t) as the magnitude (I'm not sure why) of this signal
( (1/2*sym(pi)*1j*t).*( exp(1j*1480*sym(pi)*t) - exp(1j*1440*sym(pi)*t)))
Apparently there were some missing parentheses in the denominator of the leading term, added here
( (1/(2*sym(pi)*1j*t)).*( exp(1j*1480*sym(pi)*t) - exp(1j*1440*sym(pi)*t)))
which matces h(t).
However, the ideal band pass filter has to also account for negative frequencies. So we take H(w) as defined above and add its mirror image.
H(w) = H(w) + H(-w);
Now the impulse response is
h(t) = ifourier(H(w),w,t);
[num,den] = numden(h(t));
h(t) = num/den
h(t) is basically the sum of two sinc signals, a fact that will come into play below. Also, the impulse response is non-causal, i.e., it's non-zero for t < 0, and symmetric around t = 0, facts that will come into play below. Also, it has infinite duration, i.e., it never decays exactly to zero, which will also be important.
Define x(t) and its CTFT
x(t) = 4*cos(2*sym(pi)*730*t) + 3*cos(2*sym(pi)*440*t) + 5*cos(2*sym(pi)*490*t)
X(w) = fourier(x(t),t,w)
The CTFT of the output is then
Y(w) = X(w)*H(w);
And the output is
y(t) = simplify(ifourier(Y(w),w,t))
As expected, the ideal bandpass filter exactly recovers the 730 Hz component and eliminates the other two.
If we want to compute the output of the filter via the convolution integral, it would look like this.
%syms tau real
%int(x(tau)*h(t-tau),tau,-inf,inf)
I commented it out because the computation was exceeding the Answers run time limit.
Now, if we want to try to do the analysis numerically in discrete time via the FFT, we have to deal with the fact that h(t), and x(t) for that matter, is infinite duration, and non-causal and symmetric, because FFT only applies to finite duration signals. The former means we'll have to pick a window over which h(t) is large enough to matter. The latter means we need to have our time vector be symmetric around t = 0, extending to negative time, as opposed to starting at t = 0 as in the Question. For x(t), we'll also have to pick a window as well.
Looking at h(t), it falls off as 1/(t*pi), so I'll make the h(t) window be -5 <= t <= 5, i.e., 1/(5*pi) seems pretty small. I'll arbitrarily use the same window for x(t).
Assuming the 10kHz sampling frequency, the sampled time vector is
Ts = 1e-4;
nx = -50000:50000;
nh = nx;
tvals = nx*Ts;
Evaluate x(t) and h(t) at the sampling times.
xf = matlabFunction(x(t));
hf = matlabFunction(h(t));
xval = xf(nx*Ts);
hval = hf(nh*Ts);
Recall that h(t) is really the sum of sinc functions, so we have to manually deal with the divide by 0 at t = 0 (or rewrite h(t) in terms of sinc I suppose).
limit(h(t),t,0)
hval(nh==0) = 40;
Now we can evaluate samples of y(t) using the convolution sum. We need to multiply it by Ts to have it approximate the convolution integral, because h(t) and x(t) are both finite signals (i.e., don't include impulses).
yval = conv(hval,xval)*Ts; % uses 'full' as the default
If we wanted to plot the values of y in the time domain, the corresponding time vector would be
tyval = Ts*((0:(numel(yval)-1)) + nx(1) + nh(1));
tyval covers a larger range than the time window used for h(t) and x(t), which is a consequence of convolution.
tyval([1 end])
In the frequency domain we take the fft of the input samples and multiply by Ts so the result approximates the CTFT of x(t)
X = fft(xval)*Ts; fx = (0:numel(X)-1)/numel(X)/Ts;
Do the same for the samples of the impulse response
H = fft(hval)*Ts; fh = fx;
Plot the frequency responses X and H, limited to the frequencies less than 1000 Hz for clarity
figure
plot(fx,abs(X),fh,abs(H)),xlim([0 1000]),xlabel('Hz');ylabel('Amplitude')
legend('Input','Filter')
The filter has nice, near rectangular pulse at 730 Hz. The input has the spikes at the three frequencies.
The amplitude of X at 730 Hz is ~20, whereas X(w) above has a Dirac delta at that frequency.
Recall that in the discrete-time analysis we're not working with x(t), but actually a windowed version of x(t). The CTFT of the 730 Hz component of the windowed signal is
Xr(w) = fourier(4*cos(2*sym(pi)*730*t)*rectangularPulse(-5,5,t),t,w)
And its amplitude at 730 Hz is
limit(Xr(w),2*sym(pi)*730)
as seen in the plot.
Finally, plot the frequency response of the output samples. As expected, we essentially only get the peak at 730 Hz
Y = fft(yval)*Ts; fy = (0:numel(Y)-1)/numel(Y)/Ts;
figure
plot(fy,abs(Y)),xlim([0 1000]),xlabel('Hz');ylabel('Amplitude')
The preceding analysis focused only on amplitude. Further adjustments would need to be made in the analysis if we want Y to have correct phase for comparison to the phase of the CTFT of y(t), i.e., Y(w).
