how to do implement difference equation in matlab

Hi i am stuck with this question
Write a MATLAB program to simulate the following difference equation 8y[n] - 2y[n-1] - y[n-2] = x[n] + x[n-1] for an input, x[n] = 2n u[n] and initial conditions: y[-1] = 0 and y[0] = 1
(a) Find values of x[n], the input signal and y[n], the output signal and plot these signals over the range, -1 = n = 10.
The book has told to user filter command or filtic
my code is down kindly guide me about initial conditions

6 Comments

I have done in this way
n=0:10;
a = [8 -2 -1]; % left hand side of difference equation
b = [1 1 0]; % right hand side of difference equation
x=2.^n;
k=filter(b,a,x)
plot(n,k)
Is it correct How to cater for initial conditions y[-1]=0 and y[0]=1
moonman
moonman on 14 Nov 2011
Moved: DGM on 26 Feb 2023
I have to do
Find values of x[n], the input signal and y[n], the output signal and plot these signals over the range, -1 = n = 10.
moonman
moonman on 14 Nov 2011
Moved: DGM on 26 Feb 2023
can somebody guide me
My advice:
2. Use proper punctuation mark, e.g. comma, period, question mark.
3. Read your own question again after posting, e.g. what is x[n] = 2n u[n]?
4. Update your question with comments, not answers.

Sign in to comment.

Answers (5)

I think your b is incorrect, it should be [1 1] instead. To accommodate for the initial value of y and x, you need to translate them into the corresponding filter state. filter command is an implementation of direct form II transpose, so you can use filtic to convert y and x to the state.
Here is an example, where n runs from 1 to 10. Based on you example, x[0] is 1.
n = 1:10;
a = [8 -2 -1];
b = [1 1];
yi = [1 0];
xi = 1;
zi = filtic(b,a,yi,xi)
y = filter(b,a,2.^n,zi)
BTW, I doubt if the input is really 2.^n as this becomes unbounded very quickly. Are you sure it's not 2.^(-n)?
HTH

8 Comments

moonman
moonman on 14 Nov 2011
Moved: DGM on 26 Feb 2023
Thanks Honglei this is the question of book. I hope it will be correct.U can read it again
*Write a MATLAB program to simulate the following difference equation: 8y[n] - 2y[n-1] - y[n-2] = x[n] + x[n-1] for an input, x[n] = 2n u[n] and initial conditions: y[-1] = 0 and y[0] = 1. (Hint : Try the commands ‘filter’ or ‘filtic’). a. Find the values of x[n], the input signal and y[n], the output signal and plot these signals over the range, -1 = n = 10. (All plots must be properly labeled). * So can u help me for this
Is this solution Ok
n = -1:10;
a = [8 -2 -1]; % left hand side of difference equation
b = [1 1 0]; % right hand side of difference equation
yi = [1 0];
xi = 0;
x=2*n;
zi = filtic(b,a,yi,xi)
y = filter(b,a,2*n,zi)
subplot(2,1,1)
plot(n,y)
xlabel('--------n----------->')
ylabel('-------Values------->')
title(' Plot of output signal y over the range from n=-1 to n=10')
subplot(2,1,2)
plot(n,x)
xlabel('--------n----------->')
ylabel('-------Values------->')
title(' Plot of input signal x over the range from n=-1 to n=10')
You are on the right track in general, but there are two things I want to point out:
  1. If your x is 2*n, then x(0) is 0.
  2. You already have y[0], y[-1], x[0], x[-1], so you don't need to compute them again.
Basically you should update your code to use
n = 1:10
and
xi = 0
This way, the y you get is y(1) through y(10). You can then concatenate y(0) and y(-1) to y to form the total vector. Similar things can be done for x too.
HTH
moonman
moonman on 15 Nov 2011
Moved: DGM on 26 Feb 2023
Any more help will be appreciated Can anybody else confirm me that my code is correct
You are working on a homework question. You do the best you can and if *you* cannot find an errors in your work, then you submit your answer and take your chances.
Have you considered working the results out manually and comparing them to the computed results?
Honglei Chen
Honglei Chen on 16 Nov 2011
Moved: DGM on 26 Feb 2023
I think what Walter suggests is a great way to proceed. Or you can do a program, as Fangjun suggested to output the result using a program and see if the two approaches agree. Debug and verification is just as important as programming.
BTW I think I mention in my post that you have issue with your coefficient and that is not fixed yet.
Divya K M
Divya K M on 9 May 2020
Moved: DGM on 26 Feb 2023
output of difference equation should be discrete but the output waveform is coming continuous
use stem instead of plot... you should get the correct waveform after that

Sign in to comment.

%%DSP LAb Task 4
% Difference equation implementation in matlab
%
clc
clear all
close all
% using filter function
n=-5:1:10;
index=find(n==0);
x=zeros(1,length(n));
x(index)=1;
subplot(2,2,1)
stem(n,x)
grid on
axis tight
b=[1 0];
a=[1 -2];
y=filter(b,a,x);
subplot(2,2,2)
stem(n,y,'filled','r')
grid on
axis tight
% Now without filter function
y1=zeros(1,length(n));
for i=1:length(n)
if(n(i)<0)
y1(i)=0;
end
if (n(i)>=0)
y1(i)=2*y1(i-1)+x(i);
end
end
subplot(2,2,3)
stem(n,y1,'filled','k')
grid on
axis tight

2 Comments

Is this code for given question if so means can I explain

Why do you use negative samples in the n vector? I am trying to use only positve samples and it doesnt work

Sign in to comment.

It might be a filter. But I thought all the assignment was asking you to do is to write a for-loop to generate the y series data based on the equation and the initial conditions.

1 Comment

moonman
moonman on 14 Nov 2011
Moved: DGM on 26 Feb 2023
no book has told to use filter is my code correct?? How to cater for initial conditions y[-1]=0 and y[0]=1

Sign in to comment.

% Define the symbolic variable syms z n
% Define the Z-transform of the right-hand side 2^n rhs_z = 1 / (1 - 2*z^(-1));
% Define the left-hand side: Y(z) - 2z^(-1)Y(z) + z^(-2)Y(z) lhs_z = (1 - 2*z^(-1) + z^(-2)) * sym('Y(z)');
% Set up the equation for Y(z) eq = lhs_z == rhs_z;
% Solve for Y(z) Y_z = solve(eq, 'Y(z)');
% Simplify the expression for Y(z) Y_z_simplified = simplify(Y_z);
% Perform partial fraction decomposition on Y(z) Y_z_decomp = partfrac(Y_z_simplified, z);
% Display the decomposed Y(z) disp('Decomposed Y(z):'); disp(Y_z_decomp);
% Now take the inverse Z-transform for each term y_n = iztrans(Y_z_decomp);
% Display the time-domain solution y(n) disp('The time-domain solution y(n) is:'); disp(y_n);
% Create a numerical sequence for plotting % Define the range for n (e.g., n from 0 to 10) n_values = 0:10;
% Evaluate y(n) for each n using subs (substitute n into the expression) y_values = double(subs(y_n, n, n_values));
% Plot the solution figure;
% Plot y(n) stem(n_values, y_values, 'filled', 'LineWidth', 2); title('Time-domain solution y(n)'); xlabel('n'); ylabel('y(n)'); grid on;
% Add labels to the graph for clarity text(0, y_values(1), ['y(0) = ', num2str(y_values(1))], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
% Define the symbolic variable syms z n
% Define the Z-transform of the right-hand side 2^n rhs_z = 1 / (1 - 2*z^(-1));
% Define the left-hand side: Y(z) - 2z^(-1)Y(z) + z^(-2)Y(z) lhs_z = (1 - 2*z^(-1) + z^(-2)) * sym('Y(z)');
% Set up the equation for Y(z) eq = lhs_z == rhs_z;
% Solve for Y(z) Y_z = solve(eq, 'Y(z)');
% Simplify the expression for Y(z) Y_z_simplified = simplify(Y_z);
% Perform partial fraction decomposition on Y(z) Y_z_decomp = partfrac(Y_z_simplified, z);
% Display the decomposed Y(z) disp('Decomposed Y(z):'); disp(Y_z_decomp);
% Now take the inverse Z-transform for each term y_n = iztrans(Y_z_decomp);
% Display the time-domain solution y(n) disp('The time-domain solution y(n) is:'); disp(y_n);
% Create a numerical sequence for plotting % Define the range for n (e.g., n from 0 to 10) n_values = 0:10;
% Evaluate y(n) for each n using subs (substitute n into the expression) y_values = double(subs(y_n, n, n_values));
% Plot the solution figure;
% Plot y(n) stem(n_values, y_values, 'filled', 'LineWidth', 2); title('Time-domain solution y(n)'); xlabel('n'); ylabel('y(n)'); grid on;
% Add labels to the graph for clarity text(0, y_values(1), ['y(0) = ', num2str(y_values(1))], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');

Asked:

on 14 Nov 2011

Commented:

on 5 May 2025

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!