Error using lsim(Hfinal,Y,t)

7 views (last 30 days)
Alcides Gabriel
Alcides Gabriel on 3 Oct 2022
Answered: Walter Roberson on 3 Oct 2022
The problem is that I have been "randomly" trying values on Wp1,wp2,ws1,ws2. and they execute well,they didn't give me a complex solution.
But the real cause is the lsym.I have been understanding how it work.But It seems that the problem persist.Because some type of value doesn't adress to anything.If someone can help,i will be glad
Giving me this error:
Error using DynamicSystem/lsim
When simulating the response to a specific input signal, the input data U must be a matrix with as many rows as samples in the time vector T, and as many columns as
input channels.
Error in chevisbov_I (line 49)
y=lsim(Hfinal,Y,t);
On this file:
clc; clear all; close all;
%Filtro CHEBYSHEV - RECHAZA BANDA
syms s
[Y,Fs] = audioread('senal_TP_202202_EL63 (2).wav');
wp1=2*pi*0.85*10^4;%Frecuencia de paso inferior deseada
wp2=2*pi*2.25*10^4;%Frecuencia de paso superior deseada
ws1=2*pi*1.10*10^4;%Frecuencia de rechazo inferior deseada
ws2=2*pi*1.75*10^4;%Frecuencia de rechazo superior deseada
Ap=2; %factor de atenuacion en dB de la banda de paso
As=35; %factor de atenuacion en dB de la banda de rechazo
wsprima=min(((wp2-wp1)*ws1)/(wp1*wp2-(ws1^2)) , ((wp2-wp1)*ws2)/((ws2^2)-wp1*wp2)); %ws de filtro rechaza banda
ep=sqrt((10^(Ap/10))-1); %CALCULANDO E
delta=sqrt(10^(-As/10)); %CALCULANDO DELTA
d=ep/(sqrt((delta^-2)-1)); %CALCULANDO EL PARAMETRO DISCRIMINANTE
N=3;
K=1/(cosh((acosh(d^-1))/N)); %CALCULANDO K(PARAMETRO DE SELECTIVIDAD con N=2
wpp=K*wsprima; %CALCULANDO EL WPP
N=ceil((acosh(d^-1))/(acosh(K^-1))); %VERIFICANDO LA DETERMINACION DEL ORDEN
phi=(1/N)*log((1+(1+ep^2)^(1/2))/ep); %calculando phi
c=1; %creando variable
dz=1;
for x=1:N
p(x)= (-wpp*sinh(phi)*sin(((2*x-1)*pi)/(2*N)))+(1j*wpp*cosh(phi)*cos(((2*x-1)*pi)/(2*N)));
p(x)= roundn(p(x),-6);
c=p(x)*c;
dz=(s-p(x))*dz;
end
if mod(N,2)== 0
c=(sqrt(1+(ep^2))*c); % si N es par
else
c=-c; % si N es impar
end
%h prototipo
H=c/dz ; %ft de un filtro chebychef
[numpr,denpr]=numden(H); %separando numerador con denominador
Hft=tf(sym2poly(numpr),sym2poly(denpr));%sym2poly extrae coeficientes y se crea la funcion de transferencia
Hft=minreal(Hft);
T=((wp2-wp1)*s)/(s^2+wp1*wp2) ; %filtro rechaza banda
%h final
Hfinal=subs(H,s,T); %substitucion t en s
[numf,denf]=numden(Hfinal);
Hfinal=tf(sym2poly(numf),sym2poly(denf)); %sym2poly extrae coeficientes y se crea la funcion de transferencia
Hfinalr=minreal(Hfinal);
t=0:1/Fs:13-(1/Fs); %tiempo de la señal
y=lsim(Hfinal,Y,t);
audiowrite('Chevy.wav',y,Fs)
[senal,Fs] = audioread("Chevy.wav");
p=audioplayer(senal,Fs);
play(p);
figure(1) %escalon
step(Hfinalr)
figure(2) %impulso
impulse(Hfinalr)
figure(3) %bode
bode(Hfinalr)
figure(4) %fase y magnitud
freqs(sym2poly(numf),sym2poly(denf))
figure(5) %señal de entrada y señal filtrada
plot(Y)
hold on
plot(y)
title('Frecuencia')
hold off

Answers (1)

Walter Roberson
Walter Roberson on 3 Oct 2022
y=lsim(Hfinal,Y,t);
In that call, Hfinal must be a dynamic system, and Y must be an array with as many rows as there are elements in t, and Y must have as many columns as there are input channels for the system.
Let us trace backwards to see whether that is satisfied.
Hfinal=tf(sym2poly(numf),sym2poly(denf)); %sym2poly extrae coeficientes y se crea la funcion de transferencia
sym2poly() has the property that it does not accept matrices, only vectors, so we can deduce from that that sym2poly(numf) and sym2poly(denf) are vectors rather than arrays -- and that means that you are creating a transfer function with one input and one output. So Y must be a matrix with one column, corresponding to one channel.
[Y,Fs] = audioread('senal_TP_202202_EL63 (2).wav');
Y is from audioread(). Is there proof in the code that Y will have exactly one column (one channel)? NO -- you do not validate that the number of channels is 1, so if senal_TP_202202_EL63 (2).wav contains stereo then Y will have two columns (two channels) and lsim would fail due to a mismatch between the transfer function and the data.
If you have no intention of permitting stereo, then put in an explicit test rejecting stereo files, or put in something that selects only one channel, or put in something that merges the channels down to a single channel (such as taking mean()).
t=0:1/Fs:13-(1/Fs);
That would be an appropriate time vector for 13 seconds of samples... well, except that because of accumulated round-off error, there might be one fewer entries than you might expect for 13 seconds. And there is the problem that as outside volunteers, we have no reason to expect that the input senal_TP_202202_EL63 (2).wav is exactly 13 seconds worth, precise to the sample.
I would suggest to you that you would be safer doing
t = (0:size(Y,1)-1)/Fs;
That will be exactly the right size no matter how long the input file was, and without any problems due to accumulated round-off error.

Tags

Products

Community Treasure Hunt

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

Start Hunting!