MATLAB Assignment (MATLAB SOFWARE REQUIRED) Essay Example

  • Category:
    Management
  • Document type:
    Assignment
  • Level:
    Undergraduate
  • Page:
    2
  • Words:
    1216

REPORT ON MATCHED FILTERS IMPLEMENTATION

Objective:

To gain an understanding on the use of matched filters to detect symbols in the presence of noise.

Signal S1 was generated and combined with the Gaussian noise using the code below in matlab.

% S1.m — signal generator.

% generates piecewise linear ECG signal of length L

% must post-smooth it with an N-point smoother:

% y = sgfilt(d, N, x), usually with d=0, and N=3,5,9, etc.

function x = ecg(L)

a0 = [0,1,40,1,0,-34,118,-99,0,2,21,2,0,0,0]; % template

d0 = [0,27,59,91,131,141,163,185,195,275,307,339,357,390,440];

a = a0 / max(a0);

d = round(d0 * L / d0(15)); % scale them to fit in length L

d(15)=L;

for i=1:14,

m = d(i) : d(i+1) — 1;

slope = (a(i+1) — a(i)) / (d(i+1) — d(i));

x(m+1) = a(i) + slope * (m — d(i));

When the code was implemented, the following plot was obtained.

MATLAB Assignment (MATLAB SOFWARE REQUIRED)

After showing the effect of noise, the following code was used to generate the impulse response of the optimum filter receiver and its plot is given after.

% ipulse reponse.m — ideal highpass matched filter

% s = 1, -1 = lowpass, highpass

% N = 2M+1 = filter length (odd)

% wc = cutoff frequency in [rads/sample]

function h = dlh(s, wc, N)

M = (N-1)/2;

for k = -M:M,

h(k+M+1) = (1-s) / 2 + s * wc / pi;

h(k+M+1) = s * sin(wc * k) / (pi * k);

MATLAB Assignment (MATLAB SOFWARE REQUIRED) 1

A function called mfilter.m was then developed to determine the output of the matched filter.This was done in time domain.

% mfilter.m — filtering by cascade of 2nd order sections

% [y, W] = cas(K, B, A, W, x)

% B = Kx3 numerator matrix

% A = Kx3 denominator matrix

% W = Kx3 state matrix

% x = scalar input

% y = scalar output

% based on cas.c

function [y, W] = mfilter(K, B, A, W, x)

for i = 1:K,

[y, W(i,:)] = sos(B(i,:), A(i,:), W(i,:), y);

The noisy signal was then cleaned using the following matlab code.

% sgfilt.m — filtering with length-N order-d SG smoother.

% y = sgfilt(d, N, x);

% x and y are L-dimensional column vectors; and N = 2M+1. Must have L > N+1.

% B(:, i) = b_{i-M-1} = input-on transient filters, i=1:M+1

% B(:, M+1) = b_0 = steady-state filter

% B(:, M+1+i) = b_{i} = input-off transient filters, i=0:M

function y = sgfilt(d, N, x)

M = (N-1)/2;

[L, L1] = size(x);

B = sg(d, N); % design filter

for i = 1:M+1, % input-on transients

y(i,1) = B(:,i)’ * x(1:N);

for n = M+2:L-M-1, % steady-state

y(n,1) = B(:,M+1)’ * x(n-M:n+M);

for i = 0:M, % input-off transients

y(L-M+i,1) = B(:,M+1+i)’ * x(L-N+1:L);

Part 2: Introducing synchronization

The impulse response to detect my preamble was generated by the following code;

load StudentPreamble.mat

mypreamble=Studen% parmeq.m — second-order parametric EQ filter design

% [b, a, beta] = parmeq(G0, G, GB, w0, Dw)

% b = [b0, b1, b2] = numerator coefficients

% a = [1, a1, a2] = denominator coefficients

% G0, G, GB = reference, boost/cut, and bandwidth gains

% w0, Dw = center frequency and bandwidth in [rads/sample]

% beta = design parameter

% for plain PEAK use: G0=0, G=1, GB=1/sqrt(2)

% for plain NOTCH use: G0=1, G=0, GB=1/sqrt(2)

function [b, a, beta] = parmeq(G0, G, GB, w0, Dw)

beta = tan(Dw/2) * sqrt(abs(GB^2 — G0^2)) / sqrt(abs(G^2 — GB^2));

b = [(G0 + G*beta), -2*G0*cos(w0), (G0 — G*beta)] / (1+beta);

a = [1, -2*cos(w0)/(1+beta), (1-beta)/(1+beta)];

tPreamble(76,:);

MATLAB Assignment (MATLAB SOFWARE REQUIRED) 2

The code below was used to generate a signal with a delay of 1000 samples of zeros, followed by the preamble and then symbol 1 of a thousand samples.

% sigav.m — signal averaging

% y = sigav(D, N, x)

% D = length of each period

% N = number of periods

% x = row vector of length at least ND (doesn’t check it)

% y = length-D row vector containing the averaged period

% It averages the first N blocks in x

function y = sigav(D, N, x)

for i=0:N-1,

y = y + x((i*D+1) : (i+1)*D); % accumulate i-th period

The code generated the plot below.

MATLAB Assignment (MATLAB SOFWARE REQUIRED) 3

The matched filter was then used to locate the start of the symbol. The transmitted symbol was decoded by passing it through a low pass filter and then through a reconstructor to obtain the actual symbol. The preamble is fairly long so as to minimize distortion of the signal by the channel noise.

The start of the symbol was identified and the preamble removed from the start and the end of the symbol using the following matlab code.

Part 3: ASCII Detection

load ASCIIDetection.mat

mysignal=ASCIIDetection(76,:);

function decode(File,File2)

%DECODE Converting Binary numbers to ASCII text

% DECODE(FILE,FILE2) where FILE is an input text file consisting of binary numbers

% and FILE2 is the output ASCII text file.

[spc_ent,C0]=textread(‘code.m’,’%s %s’,2);

[L,C]=textread(‘code.m’,’%s %s’,’headerlines’,2);

L=char(L);

bit = length(char(C(1)));

FID = fopen(File,’r’);

OUT = fopen(File2,’w’);

tline = fgetl(FID);

if ~ischar(tline), break, end

for i=1:size(tline,2)/bit

for j=1:93

x=isequal(char(C(j,:)),tline(8*i-7:8*i));

fprintf(OUT,’%s’,L(t));

elseif tline(8*i-7:8*i)==char(C0(2))

fprintf(OUT,’n’);

elseif tline(8*i-7:8*i)==char(C0(1))

fprintf(OUT,’ ‘);

fclose(‘all’);

The following plot of 7000 samples was obtained;

MATLAB Assignment (MATLAB SOFWARE REQUIRED) 4

The block diagram below is used to describe the detection of antipodal binary data using the matched filter.

MATLAB Assignment (MATLAB SOFWARE REQUIRED) 5

A function msampler.m was then developed to show the sampling test statistic for every symbol.

f = @(x) exp(-x.^2/2).*(1 + (sin(3*x)).^2).*…

(1 + (cos(5*x).^2));

area = quad(f,-5,5);

x = slicesample(1,N,’pdf’,f,’thin’,5,’burnin’,1000);

[binheight,bincenter] = hist(x,50);

h = bar(bincenter,binheight,’hist’);

set(h,’facecolor’,[0.8 .8 1]);

xd = get(gca,’XLim’);

xgrid = linspace(xd(1),xd(2),1000);

binwidth = (bincenter(2)-bincenter(1));

y = (N*binwidth/area) * f(xgrid);

plot(xgrid,y,’r’,’LineWidth’,2)

A function mfdecision.m was then implemented to determine the sequence of binary 1 and binary 0 transmitted.The function is given below.

noise = randn(1,size(y,2)); % random noise

ey = y + noise; % samples with noise

eY = fft(ey); % Fourier transform of noisy signal

n = size(ey,2)/2; % use size for scaling

amp_spec = abs(eY)/n; % compute amplitude spectrum

figure % plots in new window

subplot(2,1,1); % first of two plots

plot(t,ey); grid on % plot noisy signal with grid

axis([0 et -8 8]); % scale axes for viewing

xlabel(‘Time (s)’); % time expressed in seconds

ylabel(‘Amplitude’); % amplitude as function of time

subplot(2,1,2); % second of two plots

freq = (0:79)/(2*n*dt); % abscissa viewing window

plot(freq,amp_spec(1:80)); grid on % plot amplitude spectrum

xlabel(‘Frequency (Hz)’); % 1 Herz = number of cycles per second

ylabel(‘Amplitude’); % amplitude as function of frequency

Finally, a function corr_rec.m was implemented to describe a correlation receiver.

% EE 3025 Spring 2011 Discussion 4

iterations = 100000; % Lots of iterations

EXY_estimate = 0; %This is the estimate of E[XY]

EY2_estimate = 0; %This is the estimate of E[Y^2]

for iter = 1:iterations

% STEP 1: Generate gaussian random variable X with mean 0 and variance 1

X = randn;

% STEP 2: Generate gaussian random variable Z with mean 0 and variance 0.25

Z = sqrt(0.25) * randn;

% STEP 3: Compute Y from X and Z

% STEP 4: Compute X*Y, and add that to EXY_estimate

EXY_estimate = EXY_estimate + X*Y;

% STEP 5: Compute Y^2, and add that to EY2_estimate

EY2_estimate = EY2_estimate + Y^2;

%Divide by number of iterations to get the estimates

EXY_estimate = EXY_estimate/iterations;

EY2_estimate = EY2_estimate/iterations;

%Estimate C, and display it

C = EXY_estimate/EY2_estimate;

display(C);