Code

close all;
clear all;
clc;

%Filter specifications
Ap = 0.05;
Aa = 53;
OmegaP1 = 1000;
OmegaP2 = 1550;             % Index : 170286U
OmegaA1 = 1100;
OmegaA2 = 1400;
OmegaS = 3800;

TS = 2*pi / OmegaS;     % sampling period
Bt = min((OmegaA1 - OmegaP1),(OmegaP2 - OmegaA2));  %transiton bandwidth
OmegaC1 = OmegaP1 + Bt/2;   % cutoff 1
OmegaC2 = OmegaP2 - Bt/2;   % cutoff 2

deltaP = (10^(0.05*Ap) - 1)/(10^(0.05*Ap) + 1);
deltaA = 10^(-0.05*Aa);
delta = min(deltaP , deltaA);

Aa1 = -20*log10(delta);     %stop band attenuation

%calculating alpha
if (Aa1<=21)
    alpha = 0;
elseif (21<Aa1 && Aa1<=50)
    alpha = 0.5842*(Aa1 - 21)^0.4 + 0.07886*(Aa1 - 21);
else
    alpha = 0.1102*(Aa1 - 8.7);
end

%calculating D
if (Aa1<=21)
    D = 0.9222;
else
    D = (Aa1 - 7.95)/14.36;
end

%calculating N
N = ceil(OmegaS*D/Bt + 1);
if(mod(N,2)==0)
    N = N+1;
end


%calculations of window function
halfRange = (N-1)/2;
n = -halfRange : 1 : halfRange;
beta = alpha*(1 - (2*n/(N-1)).^2).^0.5;

%Generating Io_alpha
bessellimit = 125;
Io_alpha = 1;
for k = 1:bessellimit
    val_k = ((1/factorial(k))*(alpha/2).^k).^2;
    Io_alpha = Io_alpha + val_k;
end

%Generating Io_beta
Io_beta = 1;
for m = 1:bessellimit
    val_m = ((1/factorial(m))*(beta/2).^m).^2;
    Io_beta = Io_beta +val_m;
end

wk_nT = Io_beta/Io_alpha;
%Printing the results
% fprintf('Filter order = %d',N);

%Plotting the kaiser function
figure;
stem(n,wk_nT);
xlabel('n');
ylabel('Amplitude');
title('Kaiser window(Time domain)');

%h[n]
neg = -halfRange : 1 : -1;
hneg = ((1/pi)./neg).*(sin(OmegaC1*TS.*neg) - sin(OmegaC2*TS.*neg));
hzero = 1 + 2*(OmegaC1 - OmegaC2)/OmegaS;
nposi = 1 : 1 : halfRange;
hposi = ((1/pi)./nposi).*(sin(OmegaC1*TS.*nposi) - sin(OmegaC2*TS.*nposi));
h = [hneg,hzero,hposi];
n = [neg,0,nposi];

figure;
stem(n,h);
grid on;
xlabel('n');
ylabel('h[n]');
title('Impulse Response');

%Filter response
hw_nT = h.*wk_nT;

%Question 2
%Plotting the causal stopband filter
n_shifted = [0:1:N-1];
figure;
stem(n_shifted,hw_nT);
xlabel('n');
ylabel('Amplitude');
title('Causal Impulse Response filter(Time Domain)');

%obtaining the frequency domain impulse response response
fvtool(hw_nT);

%Question 3
[Hw,f] = freqz(hw_nT);      %obtaining the frequency response and corresponding frequencies
w = f*OmegaS/(2*pi);        %Angular frequency
log_Hw = 20.*log10(abs(Hw));
figure;
plot(w,log_Hw);
xlabel('Angular frequency(rad/s)');
ylabel('Magnitude(dB)');
title('Magnitude response of the filter(Frequency domain)');

%Question 4
%Plotting the magnitude response of the lower passband
figure;
finish = round((length(w)/(OmegaS/2)*OmegaC1));
wpass_l = w(1:finish);
hpass_l = log_Hw(1:finish);
plot(wpass_l,hpass_l);
axis([-inf, inf, -0.1, 0.1]);
xlabel('Frequency (rad/s)');
ylabel('Magnitude (dB)');
title('Magnitude response of Lower Passband - Frequency Domain');

%Plotting the magnitude response of the lower passband
figure;
start = round(length(w)/(OmegaS/2)*OmegaC2);
wpass_h = w(start:length(w));
hpass_h = log_Hw(start:length(w));
plot(wpass_h,hpass_h);
axis([-inf, inf, -0.1, 0.1]);
xlabel('Frequency (rad/s)');
ylabel('Magnitude (dB)');
title('Magnitude response of the Upper Passband - Frequency Domain');


%Question 5
%Generating the input of desired number samples
%Generating the discrete signal
%Component frequencies of the input
O_1 = OmegaC1/2;
O_2 = OmegaC1 + (OmegaC2-OmegaC1)/2;
O_3 = OmegaC2 + (OmegaS/2-OmegaC2)/2;

n1 = 0:1:600;   % samples
X = cos(O_1.*n1.*TS)+cos(O_2.*n1.*TS)+cos(O_3.*n1.*TS);

figure;
subplot(2,1,1);
stem(n1,X);
xlabel('n');
ylabel('Amplitude');
title('Input signal(Time domain)')
subplot(2,1,2);
len_fft = 2^nextpow2(numel(n1))-1;
x_fft = fft(X,len_fft);
x_fft_plot = [abs([x_fft(len_fft/2+1:len_fft)]),abs(x_fft(1)),abs(x_fft(2:len_fft/2+1))];
f = OmegaS*linspace(0,1,len_fft)-OmegaS/2;
plot(f,x_fft_plot);
xlabel('Frequency rad/s');
ylabel('Magnitude');
title('Input signal in the frequency domain');
axis tight;

%Question 6
% Filtering using frequency domain multiplication
len_fft = length(X)+length(hw_nT)-1; % length for fft in x dimension
x_fft = fft(X,len_fft);
hw_nT_fft = fft(hw_nT,len_fft);
out_fft = hw_nT_fft.*x_fft;
out = ifft(out_fft,len_fft);
rec_out = out(floor(N/2)+1:length(out)-floor(N/2));

% Ideal Output Signal
ideal_out = cos(O_1.*n1.*TS)+cos(O_3.*n1.*TS);

%Obtaining the output waveforms
% Frequency domain representation of output signal after filtering using the designed filter
figure;
subplot(2,1,1);
len_fft = 2^nextpow2(numel(n1))-1;
xfft_out = fft(rec_out,len_fft);
x_fft_out_plot = [abs([xfft_out(len_fft/2+1:len_fft)]),abs(xfft_out(1)),abs(xfft_out(2:len_fft/2+1))];
f = OmegaS*linspace(0,1,len_fft)-OmegaS/2;
plot(f,x_fft_out_plot);
xlabel('Frequency rad/s');
ylabel('Magnitude');
title('Output signal of the designed filter in the frequency domain');

% Time domain representation of output signal after filtering using the designed filter
subplot(2,1,2);
stem(n1,rec_out);
xlabel('n');
ylabel('Amplitude');
title('Output signal of the designed filter in the time domain');

%Obtaining the outputs of the ideal filter
figure;
subplot(2,1,1);
xfft_outideal = fft(ideal_out,len_fft);
x_fft_outideal_plot = [abs([xfft_outideal(len_fft/2+1:len_fft)]),abs(xfft_outideal(1)),abs(xfft_outideal(2:len_fft/2+1))];
plot(f,x_fft_outideal_plot);
xlabel('Frequency rad/s');
ylabel('Magnitude');
title('Output signal of the ideal filter in the frequency domain');

% Time domain representation of output signal after filtering using ideal filter
subplot(2,1,2);
stem(n1,ideal_out);
xlabel('n');
ylabel('Amplitude');
title('Output signal of the ideal filter in the time domain');