Lab Experiment-1 (A) Conversion of A Matrix Into Its Row Echelon Form
Lab Experiment-1 (A) Conversion of A Matrix Into Its Row Echelon Form
Lab Experiment-1 (A) Conversion of A Matrix Into Its Row Echelon Form
Aim: To convert a matrix into its row echelon form Theory: We can convert any matrix into its row reduced echelon form by applying elementary row operations. Elementary row operations preserve the row space of the matrix. Procedure for getting the Row echelon form of a matrix A matrix is said to be in row echelon form if the following conditions are satisfied. 1) All nonzero rows (rows with at least one nonzero element) are above any rows of all zeroes (All zero rows, if any, belong at the bottom of the matrix). 2) The leading coefficient (the first nonzero number from the left, also called the pivot) of a nonzero row is always strictly to the right of the leading coefficient of the row above it. 3) All entries in a column below a leading entry are zeroes.
Matlab code clear all; close all; clc; a =input('enter the matrix'); [r c]=size(a); for j=1:r; for i=1:r-1; if (a(i,1)<a(i+1,1)) a([i,i+1],:)=a([i+1,i],:);%sorting of rows end end end for i=1:min(r,c) if (a(i,i)~=0) a(i,:)=(1/a(i,i)).*a(i,:);%normalising the pivot rows end end for j=1:min(c,r) for i=1:min(c,r)
if(i~=j) a(i,:)=a(i,:)*a(j,j)-a(j,:)*a(i,j);%echolon form end end end for i=1:min(r,c) if (a(i,i)~=0) a(i,:)=(1/a(i,i)).*a(i,:);%normalising the pivot rows end end Program Output 1) Using the designed code a= 1.0000 0 0 -0.2222 0 1.0000 0 3.0889 0 0 1.0000 -1.9556 2) Using the built in Matlab command (rref) rref([1 2 1 4;2 3 4 1;10 2 1 2]) ans = 1.0000 0 0 -0.2222 0 1.0000 0 3.0889 0 0 1.0000 -1.955
Lab Experiment-1(b) To implement Interpolation and Decimation functions Aim: Set up a program to up sample (Interpolation) and down sample (Decimation) the signal by a given factor. Theory: Interpolation Interpolation is also known as up sampling. It is the process of increasing the sampling rate of the given signal. Interpolation causes signal expansion in time domain where as signal compression in frequency domain. This process is mainly used in multirate signal processing. The up sampling factor (commonly denoted by L) is usually an integer or a rational fraction greater than unity. This factor multiplies the sampling rate or, equivalently, divides the sampling period. Interpolation is defined by the function:
n 0, L, 2L, otherwise
Decimation Decimation or down sampling means reducing the sampling rate. It implies compression in time domain and expansion in frequency domain. The down sampling factor (commonly denoted by M) is usually an integer or a rational fraction greater than unity. Down sampling operation is implemented by keeping every Mth sample and removing all M-1 samples to generate xint[n]. Down sampling is denoted by the equation:
Matlab Code (For Interpolator) y=ones(1,8); L=input( 'enter the interpolation factor') x=zeros(1,length(y)*L); x(1:L:length(x))=y; N=50;%order of the filter n=0:N-1; M=N/2; wc=pi/L;%cutoff frequency for i=1:length(n) if (n(i)~=M) Hlp(i)=(sin(wc*(n(i)-M)))/(pi*(n(i)-M)); else Hlp(i)=wc/pi; end end %creating Hamming Window for i=1:length(n) w(i)=0.54+0.46*(1-cos((2*pi*n(i))/(N-1))); end %multiplying the ideal response with the window % %---------------------------------------------------------% %desired filter response Hd=Hlp.*w; yout=conv(Hd,x) wn = -1.5*pi:0.01:1.5*pi; X=freqz(y,1,wn); Y=freqz(x,1,wn); Z=freqz(yout,1,wn); subplot(3,1,1) plot(wn,abs(X)) xlabel('w') ylabel('magnitude') title('Original signal') subplot(3,1,2) plot(wn,abs(Y)) xlabel('w') ylabel('magnitude') title('O/p of interpolator') subplot(3,1,3) plot(wn,abs(Z)) xlabel('w') ylabel('magnitude') title('Filtered interpolator output')
figure subplot(3,1,1) stem(y) xlabel('n') ylabel('magnitude') title('Original signal') subplot(3,1,2) stem(x) title('O/p of interpolator') xlabel('n') ylabel('magnitude') subplot(3,1,3) stem(yout) xlabel('n') ylabel('magnitude') title('Filtered interpolator output')
Results
Original signal
magnitude
1 0.5 0 1 2 3 4 5 6 7 8
n O/p of interpolator
magnitude
magnitude
1 0 -1 0 10 20 30 40 n 50 60 70 80
Original signal
magnitude
10 5 0 -5 -4 -3 -2 0 1 w O/p of interpolator -1 2 3 4 5
magnitude
magnitude
20 10 0 -5 -4 -3 -2 -1 0 w 1 2 3 4 5
Figure-2: Frequency domain plots Matlab Code (For Decimator) clc close all x=ones(1,8); n=0:length(x)-1; M=input('enter the decimation factor:') y=[ ]; k=[ ]; for i1=1:length(n) n1=n(i1); if (mod(n1,M)==0) k=[k n1]; y=[y x(i1)]; end end k=(1/M).*k;
wn=-1.5*pi:0.01:1.5*pi; X=freqz(x,1,wn); Y=freqz(y,1,wn); subplot(2,1,1) plot(wn,abs(X)) xlabel('w') ylabel('magnitude') title('Original signal') subplot(2,1,2) plot(wn,abs(Y)) xlabel('w') ylabel('magnitude') title('O/p of decimator') figure subplot(2,1,1) stem(n,x) xlabel('n') ylabel('x[n]') legend('Input signal') subplot(2,1,2) stem(k,y) xlabel('k') ylabel('y[k]') legend('Output signal')
Results
1 Input signal
x[n]
0.5
3 n
1 Output signal
y[k]
0.5
0.5
1.5 k
2.5
Original signal 8
magnitude
6 4 2 0 -5
-4
-3
-2
-1
0 1 w O/p of decimator
magnitude
3 2 1 0 -5
-4
-3
-2
-1
0 w
Discussion
From the above figures ,we can clearly see the time and frequency domain plots of interolation and decimation and the phenomenon of Imaging in Interpolation.It has been eliminated through use of an suitable anti imaging filter.
Lab Experiment-2 Design of FIR digital filters Aim: To design Low Pass, High Pass, Band Pass and Band Stop FIR digital filters Theory: The impulse responses of ideal low-pass, high-pass, band-pass and band-stop filters are as given in the following table.
Since these impulse responses are of infinite duration, we multiply it with a suitable window function to limit it to a finite length. The various characteristics of the filter like transition region width, pass-band and stop-band attenuation etc depend upon the window function which is selected. In the following MATLAB programs, we have used the Hann window, Hamming window, and Blackmann window for the implementation of the different types of filters mentioned above, according to user defined specifications.
1) Low-pass filter Matlab code clc rp=input('enter passband ripple(in dB):'); rs=input('enter stop band ripple (in dB):'); fp=input('passband frequency:'); fs=input('stopband frequency:'); f=input('sampling frequency:'); rp=10^(rp/20); rs=10^(rs/20); wp=2*fp/f; ws=2*fs/f; fc=(fs+fp)/2; wcr=2*pi*fc/f; num=-20*log10(sqrt(rp*rs))-13; dem=14.6*(fs-fp)/f; n=ceil(num/dem); n1=n+1; if(rem(n,2)~=0) n1=n; n=n-1; end k=0:n1-1; M=n/2; % %---------------------------------------------------------%Calculation of ideal filter impulse response for i=1:length(k) if (k(i)~=M) Hlp(i)=(sin(wcr*(k(i)-M)))/(pi*(k(i)-M)); else Hlp(i)=(wcr)/pi; end end % %---------------------------------------------------------% %Creating Hann window for i=1:length(k) if rs<44 w(i)=0.5-0.5*cos(2*pi*k(i)/(n-1)); % Hanning window else if rs>=44 && rs<53 w(i)=0.54-0.46*cos(2*pi*k(i)/(n-1)); % Hamming window else if rs>=53
w(i)=0.42-0.5*cos(2*pi*k(i)/(n-1))+0.08*cos(4*pi*k(i)/(n-1));% Blackmann window end end end end % %---------------------------------------------------------% %desired filter response Hd=Hlp.*w; fvtool(Hd)
Results
enter passband ripple(in dB):-60 enter stop band ripple (in dB):-60 passband frequency:200/6 stopband frequency:50 sampling frequency:200
Magnitude Response (dB) 0
-10
-20
-30
Magnitude (dB)
-40
-50
-60
-70
-80
-90 0 0.1 0.2 0.3 0.4 0.5 0.6 Normalized Frequency ( rad/sample) 0.7 0.8 0.9
Phase Response 0
-10
-20
-30
Phase (radians)
-40
-50
-60
-70
-80 0 0.1 0.2 0.3 0.4 0.5 0.6 Normalized Frequency ( rad/sample) 0.7 0.8 0.9
2) Band-pass filter Matlab code clc rp=input('enter passband ripple(in dB):'); rs=input('enter stop band ripple (in dB):'); fp1=input('enter passband frequency 1:'); fs1=input('enter stopband frequency 1:'); fp2=input('enter passband frequency 2:'); fs2=input('enter stopband frequency 2:'); f=input('sampling frequency:'); rp=10^(rp/20); rs=10^(rs/20); fc1=(fs1+fp1)/2; fc2=(fs2+fp2)/2; wcr1=2*pi*fc1/f; wcr2=2*pi*fc2/f; num=-20*log10(sqrt(rp*rs))-13; dem=14.6*(fs-fp)/f; n=ceil(num/dem); n1=n+1; if(rem(n,2)~=0) n1=n; n=n-1; end k=0:n1-1; M=n/2; % %---------------------------------------------------------%Calculation of ideal filter impulse response for i=1:length(k) if (k(i)~=M) Hlp(i)=((sin(wcr2*(k(i)-M)))/(pi*(k(i)-M)))-((sin(wcr1*(k(i)-M)))/(pi*(k(i)-M))); else Hlp(i)=(wcr2-wcr1)/pi; end end % %---------------------------------------------------------% %Creating Hann window for i=1:length(k) if rs<44 w(i)=0.5-0.5*cos(2*pi*k(i)/(n-1)); % Hanning window else if rs>=44 && rs<53 w(i)=0.54-0.46*cos(2*pi*k(i)/(n-1)); % Hamming window else if rs>=53 w(i)=0.42-0.5*cos(2*pi*k(i)/(n-1))+0.08*cos(4*pi*k(i)/(n-1));% Blackmann window end
Results
enter passband ripple(in dB):-60 enter stop band ripple (in dB):-60 enter passband frequency 1:55 enter stopband frequency 1:45 enter passband frequency 2:145 enter stopband frequency 2:155 sampling frequency:400
Magnitude Response (dB) 0
-10
-20
-30
Magnitude (dB)
-40
-50
-60
-70
-80 0 0.1 0.2 0.3 0.4 0.5 0.6 Normalized Frequency ( rad/sample) 0.7 0.8 0.9
Phase Response
-20
-40
Phase (radians)
-60
-80
-100
-120 0 0.1 0.2 0.3 0.4 0.5 0.6 Normalized Frequency ( rad/sample) 0.7 0.8 0.9
3) Band-stop filter Matlab code clc rp=input('enter passband ripple(in dB):'); rs=input('enter stop band ripple (in dB):'); fp1=input('enter passband frequency 1:'); fs1=input('enter stopband frequency 1:'); fp2=input('enter passband frequency 2:'); fs2=input('enter stopband frequency 2:'); f=input('sampling frequency:'); rp=10^(rp/20); rs=10^(rs/20); fc1=(fs1+fp1)/2; fc2=(fs2+fp2)/2; wcr1=2*pi*fc1/f; wcr2=2*pi*fc2/f; num=-20*log10(sqrt(rp*rs))-13; dem=14.6*(fs-fp)/f; n=ceil(num/dem); n1=n+1; if(rem(n,2)~=0) n1=n; n=n-1; end k=0:n1-1; M=n/2; % %----------------------------------------------------------
%Calculation of ideal filter impulse response for i=1:length(k) if (k(i)~=M) Hlp(i)=((sin(wcr1*(k(i)-M)))/(pi*(k(i)-M)))-((sin(wcr2*(k(i)-M)))/(pi*(k(i)-M))); else Hlp(i)=1-((wcr2-wcr1)/pi); end end % %---------------------------------------------------------% %Creating Hann window for i=1:length(k) if rs<44 w(i)=0.5-0.5*cos(2*pi*k(i)/(n-1)); % Hanning window else if rs>=44 && rs<53 w(i)=0.54-0.46*cos(2*pi*k(i)/(n-1)); % Hamming window else if rs>=53 w(i)=0.42-0.5*cos(2*pi*k(i)/(n-1))+0.08*cos(4*pi*k(i)/(n-1));% Blackmann window end end end end % %---------------------------------------------------------% %desired filter response Hd=Hlp.*w; fvtool(Hd)
Results
enter passband ripple(in dB):-60 enter stop band ripple (in dB):-60 enter passband frequency 1:70 enter stopband frequency 1:80 enter passband frequency 2:130 enter stopband frequency 2:120 sampling frequency:400
-10
-20
Magnitude (dB)
-30
-40
-50
-60
0.1
0.2
0.3
0.7
0.8
0.9
Phase Response 0
-20
-40
Phase (radians)
-60
-80
-100
0.1
0.2
0.3
0.7
0.8
0.9
4) High-pass Filter clc rp=input('enter passband ripple(in dB):'); rs=input('enter stop band ripple (in dB):'); fp=input('enter passband frequency (in Hz):'); fs=input('enter stopband frequency (in Hz):'); f=input('enter sampling frequency (in Hz):'); rp=10^(rp/20); rs=10^(rs/20);
wp=2*fp/f; ws=2*fs/f; fc=(fs+fp)/2; wcr=2*pi*fc/f; num=-20*log10(sqrt(rp*rs))-13; dem=14.6*(fs-fp)/f; n=ceil(num/dem); n1=n+1; if(rem(n,2)~=0) n1=n; n=n-1; end k=0:n1-1; M=n/2; % %---------------------------------------------------------%Calculation of ideal filter impulse response for i=1:length(k) if (k(i)~=M) Hlp(i)=-((sin(wcr*(k(i)-M)))/(pi*(k(i)-M))); else Hlp(i)=1-(wcr/pi); end end % %---------------------------------------------------------% %Creating Hann window for i=1:length(k) if rs<44 w(i)=0.5-0.5*cos(2*pi*k(i)/(n-1)); % Hanning window else if rs>=44 && rs<53 w(i)=0.54-0.46*cos(2*pi*k(i)/(n-1)); % Hamming window else if rs>=53 w(i)=0.42-0.5*cos(2*pi*k(i)/(n1))+0.08*cos(4*pi*k(i)/(n-1));% Blackmann window end end end end % %---------------------------------------------------------% %desired filter response Hd=Hlp.*w; fvtool(Hd)
Results
enter passband ripple(in dB):-60 enter stop band ripple (in dB):-60 enter passband frequency (in Hz):95 enter stopband frequency (in Hz):105
-40
Magnitude (dB)
-60
-80
-100
-120 0 0.1 0.2 0.3 0.4 0.5 0.6 Normalized Frequency ( rad/sample) 0.7 0.8 0.9
Phase Response
-20
-40
-60
Phase (radians)
-80
-100
-120
-140
-160
0.1
0.2
0.3
0.7
0.8
0.9
Discussion
We have sucessfully implemented Low-pass, High-pass, Band-pass and Band-stop filters using Windowing method with user defined specifications. The parameters like stop-band ripple, passband ripple, stop band and pass band frequencies are taken as input from the user. The program decides the window function based on the users inputs and designs the filters with required characteristics.
Aim: To find the best approximation to a sine wave by using polynomials of degree less than or equal to 5. Theory: Curve fitting is the process of constructing a curve or mathematical functions that has a best fit to a series of data-points. curve fitting can be used as an aid for data visualization, to infer values of a function where no data are available and also to summarize the relationships among two or more variables. For a given fifth degree polynomial, the set {1, X, X2, X3, X4, X5} forms a basis. In order to produce corresponding standard basis (e1, e2, e3, e4, e5, e6), we apply Gram-Schmidt procedure. For this vector space, the inner product is defined by the following equation: <f, g>= ( ) ( )dx
Now, In order to get the approximation, sin(x) is projected on to the standard basis of P 5 to find out the coefficients corresponding to each basis function and using them to evaluate the following expression: Pv=<v, e1>e1+<v, e2>e2+<v, e3>+<v, e4>e4+<v, e5>e5+<v, e6>e6 where Pv is the polynomial, which is the best approximation to v(x)=sin(x), v C (, ) such that the error in approximation( i.e., ||v-Pv||) is minimized.
Matlab Code % to find the best fitting polynomial of at most fifth degree clc clear %defininng set of initial basis y=sym('sin(x)'); u1=sym('1'); u2=sym('x'); u3=sym('x^2'); u4=sym('x^3'); u5=sym('x^4'); u6=sym('x^5'); %doing gram schmidt v1=u1; v2=u2-((inner_pdt(u2,v1))/(inner_pdt(v1,v1)))*v1; v3=u3((((inner_pdt(u3,v1))/(inner_pdt(v1,v1)))*v1)+((inner_pdt(u3,v2))/(inner_pdt(v2,v2)))*v2); v4=u4(((((inner_pdt(u4,v1))/(inner_pdt(v1,v1)))*v1+((inner_pdt(u4,v2))/(inner_pdt(v2,v2)))*v2+((inn er_pdt(u4,v3))/(inner_pdt(v3,v3)))*v3))); v5=u5(((inner_pdt(u5,v1))/(inner_pdt(v1,v1)))*v1+((inner_pdt(u5,v2))/(inner_pdt(v2,v2)))*v2+((inner _pdt(u5,v3))/(inner_pdt(v3,v3)))*v3+((inner_pdt(u5,v4))/(inner_pdt(v4,v4)))*v4); v6=u6(((inner_pdt(u6,v1))/(inner_pdt(v1,v1)))*v1+((inner_pdt(u6,v2))/(inner_pdt(v2,v2)))*v2+((inner _pdt(u6,v3))/(inner_pdt(v3,v3)))*v3+((inner_pdt(u6,v4))/(inner_pdt(v4,v4)))*v4+((inner_pdt(u6 ,v5))/(inner_pdt(v5,v5)))*v5); v11=v1/(inner_pdt(v1,v1)^0.5); v22=v2/(inner_pdt(v2,v2)^0.5); v33=v3/(inner_pdt(v3,v3)^0.5); v44=v4/(inner_pdt(v4,v4)^0.5); v55=v5/(inner_pdt(v5,v5)^0.5); v66=v6/(inner_pdt(v6,v6)^0.5); Pu_v=(inner_pdt(y,v11)*v11)+(inner_pdt(y,v22)*v22)+(inner_pdt(y,v33)*v33)+(inner_pdt(y,v4 4)*v44)+(inner_pdt(y,v55)*v55)+(inner_pdt(y,v66)*v66) ezplot(Pu_v,[-pi,pi]) figure ezplot(y,[-pi,pi]) inner_pdt.m
function [y]=inner_pdt(x1,x2) fun=sym(x1*x2); y=int(fun,-pi,pi);
Results
Plot of Sin(x) vs x
0.5
Sin(x)
-0.5
-1
-3
-2
-1
0 x
0.5
Puv
-0.5
-1
-3
-2
-1
0
x
Figure-2: Plot of the best approximation to Sin(x) using 5th degree Polynomial
Figure-2: Typical frequency response plots of the filters used in a 2 channel QMF bank
Matlab Code % To implement the quadrature mirror filter bank clc clear %designing the input signal f=1; t=0:1/100:1; x_n=sin(2*pi*f*t); %-------------------------------------------------------------------------% design of Ho(Z)-Prototype filter disp('Enter the specifications of prototype filter Ho(z)') rp=input('Enter passband ripple(in dB):'); rs=input('Enter stop band ripple (in dB):'); fp=input('Enter passband frequency (in Hz):'); fs=input('Enter stopband frequency (in Hz):'); f=input('Enter sampling frequency (in Hz):'); rp=10^(rp/20); rs=10^(rs/20); wp=2*fp/f; ws=2*fs/f; fc=(fs+fp)/2; wcr=2*pi*fc/f; num=-20*log10(sqrt(rp*rs))-13; dem=14.6*(fs-fp)/f; N=ceil(num/dem); n1=N+1; if(rem(N,2)~=0) n1=N; N=N-1; end n=0:N-1; M=N/2; % %---------------------------------------------------------%Calculation of ideal filter impulse response for i=1:length(n) if (n(i)~=M) Hlp(i)=(sin(wcr*(n(i)-M)))/(pi*(n(i)-M)); else Hlp(i)=(wcr)/pi; end end % %---------------------------------------------------------% %Creating Hann window for i=1:length(n) if rs<44
w(i)=0.5-0.5*cos(2*pi*n(i)/(N-1)); % Hanning window else if rs>=44 && rs<53 w(i)=0.54-0.46*cos(2*pi*n(i)/(N-1)); % Hamming window else if rs>=53 w(i)=0.42-0.5*cos(2*pi*n(i)/(N-1))+0.08*cos(4*pi*n(i)/(N-1));% Blackmann window end end end end % %desired filter response Ho=Hlp.*w; %-------------------------------------------------------------------------%-------------------------------------------------------------------------% designing H1(Z) cmplx_expo=exp(-j*pi.*n); H1=Ho.*cmplx_expo; %-------------------------------------------------------------------------%top branch of the system x2=conv(x_n,Ho); x3=decimate(x2,2); x4=interp(x3,2); x5=conv(x4,Ho); %-------------------------------------------------------------------------%lower branch of the system v2=conv(x_n,H1); v3=decimate(v2,2); v4=interp(v3,2); v5=conv(v4,-1.*H1); %-------------------------------------------------------------------------%adding the output of both the branches y1=v5+x5; % compensation for delay introduced by Filters n1=1:length(y1); [y2 k2]=sig_shift(y1,n1,-2*M); y=y2((2*M+1):(2*M+1)+(length(x_n)-1)); k=k2((2*M+1):(2*M+1)+(length(x_n)-1)); subplot(3,1,1) stem(x2) ylabel('Vo[n]') xlabel('n') title('Output of top lowpass filter (Ho(Z))') subplot(3,1,2) stem(x3) ylabel('Uo[n]')
xlabel('n') title('Output of top Decimator') subplot(3,1,3) stem(x3) ylabel('V^o[n]') xlabel('n') title('Output of top Interpolator') figure subplot(3,1,1) stem(v2) ylabel('V1[n]') xlabel('n') title('Output of highpass filter (H1(Z))') subplot(3,1,2) stem(v3) ylabel('U1[n]') xlabel('n') title('Output of bottom Decimator') subplot(3,1,3) stem(v4) ylabel('V^1[n]') xlabel('n') title('Output of bottom Interpolator') figure subplot(2,1,1) stem(x_n) legend('input signal') ylabel('x[n]') xlabel('n') subplot(2,1,2) stem(k,y) legend('output signal') ylabel('y[n]') xlabel('k')
Results
Enter the specifications of prototype filter Ho(z) Enter passband ripple(in dB):-60 Enter stop band ripple (in dB):-60 Enter passband frequency (in Hz):95 Enter stopband frequency (in Hz):105
Vo[n]
Uo[n]
Vo[n]
x 10
-15
V1[n]
0 -2 0 x 10
-15
50
100
200
250
U1[n]
0 -2 0 x 10
-15
20
40
100
120
V1[n]
x[n]
0 -0.5 -1
20
40
60 n
80
100
120
y[k]
0 -0.5 -1
20
40
60 k
80
100
120
Figure-3: Input Signal and Shifted version of Output Signal Discussion The output signal is only an amplitude-scaled and time shifted version of the input signal. We have compensated for this time shift by adding appropriate time advance blocks just for the purpose of clarity in the plotting output signal. Thus, Here we have obtained Perfect Reconstruction.