%HELSINKI UNIVERSITY OF TEHNOLOGY
%Communications Laboratory
%Seppo Saastamoinen 5.2.2005

%S-72.060 SIGNALS AND SYSTEMS
%Laboratory work 2: Examing the filters
%Makro: filter1.m

clear all; close all; %Clear all variables and figures.

%With this macro you can design filters.

%Basic values:

Wp = [1000 3000]; %Band pass corner frequencies.
Ws = [100 3900]; %Stop band corner frequencies.

[n,Wn] = buttord(Wp,Ws,3,20,'s'); %Filter's grade n and normalized
%pass band frequencies.

%Below is calculated the transfer function´s H(s) numerator b(s) and denominator a(s)
%using Laplace transform where s = jw.

[b,a] = butter(n,Wn,'s'); %Transfer function H(s)=b(s)/a(s).

f = (0:1:20000); %Setting the frequency domain.

%Below the frequency response and angle function are calculated:
Hf = freqs(b,a,f);

%Below the amplitude response in dB´s is calculated:
Af = 20*log10(abs(Hf)); %Af is the amplification of the filter.

%Below the radians are changed to degrees and angle function is on the area
%-180...180 degrees:
Pf = (180/pi)*angle(Hf);

%Below the angle accumulation is calculated:
P1 = unwrap(angle(Hf));

%Below the differential coefficient of frequency is calculated:
df = diff(f);

%Below the differential coefficient of angle accumulation is calculated:
dP = diff(P1);

%Below the group delay is calculated:
G = (-1/(2*pi))*dP./df;

%Below the lenght of the vector G is corrected to have same the lenght as other vectors:
Gd = [0 G];

%Setting the plotting section:.
subplot(3,1,1); %3 figures.

semilogx(f,Af); %Plotting the amplitude function (response).
%Y axis is linear, X axis is logarithmic.
%If needed X axis can be linear. The use command plot instedt of semilogx.

%The title for the page can have many rows as below::
title ({'S-72.060 Signals and ja järjestelmät, laboratory works 2005',...
      'Work 2: Examing the filters','Butterworth band pass filter with N = 3, 1000-3000 Hz',...
      'Assistant: N.N','Students: Niina Näppärä 65543G ja Oiva Piirakainen 23777K 5.2.2005'});

grid on; %Seting the grid on.
axis([100 20000 -120 3]); %Setting the axis limits.
xlabel('Frequency/Hz'); %Name of the x axis.
ylabel('Gain/dB'); %Name of the y axis.
zoom on; %Zooming use mouse. Double click returns the original figure.

subplot(3,1,2);
plot(f,Pf); %Plotting the angle response.
grid on;
axis([100 20000 -200 200]);
xlabel('Frequency/Hz');
ylabel('Angle/degrees');
zoom on;

subplot(3,1,3);
plot(f,Gd); %Plotting the group delay.
grid on;
axis([100 20000 -0.0001 0.005]);
xlabel('Frequency/Hz');
ylabel('Group delay/s');
zoom on;

shg; %Allows the figure to be at the top of the screen.
