Добрый вечер, уважаемые форумчане! Прошу помощи, необходимо в системе Matlab повторить все аналогичные действия, перечисленные в методичке пункты 2-5.
Проблема в том, что не могу построить обратное преобразование Фурье и вывести его, а также построить корреляционные (АКФ) функции сигналов. Также при записи в согласованной фильтрации, при ее выполнении для экспоненциально затухающего импульса и прямоугольного импульса, начинающихся в середине интервала, на выходе фильтра получается не то что надо (для экспоненциально затухающего импульса на выходе должна быть линия, возрастающая после начала импульса в середине, а для прямоугольного импульса должен получиться график в виде треугольника, максимум которого достигается в точке начала следования прямоугольного импульса).
Я повторил аналогичные действия в Mathcad, чтобы визуализировать то, что надо сделать в матлабе. На всякий случай скину методичку и сами скриншоты из маткада. Не получается повторить все те же действия, что и в маткаде, но с помощью матлаб... Не могу сделать, даже за деньги заказ не берут, а осталось всего-то, но вот не могу исправить ошибки, которые выдает программа! Буду благодарен за оказанную помощь!
Методические указания :
https://drive.google.com/file/d/16MEhu72rEczrhVsHcTUV5OXZCMST1EN-/view?usp=sharingВот код в матлаб:
% % узкополосная фильтрация
%
x0 = 1.5;
w0 = 2*pi*60;
t = [0:0.001:0.8]';
N = 4096;
n = length(t);
Y = zeros(1,n);
for i = 1:n
if t(i) >= 0 && t(i) <= 0.4
Y(i) = 0;
elseif t(i) > 0.4 && t(i) <= 0.8
Y(i) = x0*cos(w0*t(i));
end
end
plot (t,Y)
% проводим дискретизацию по времени
n1 = [0:N-1];
tt = n1.*(t(end)/N);
for i = 1:length(n1)
if tt(i) >= 0 && tt(i) <= 0.4
Y1(i) = 0;
elseif tt(i) > 0.4 && tt(i) <= 0.8
Y1(i) = x0*cos(w0*tt(i));
end
end
figure,plot(tt,Y1)
% накладываем гауссовский белый шум
sigma = 2;
Z = Y1+normrnd(0,sigma,size(n1));
figure,plot(tt,Z,'k-')
% подвергаем преобразованию Фурье
%L = [1:N-1]';
%nfft = length(L);
ns = (length(tt));
ks=1/tt(2);
YY = (abs(fft(Y1))/(ns));
YY=YY(1:ns/2+1);
YY(2:end-1) = 2*YY(2:end-1);
ZZ = (abs(fft(Z))/(ns));
ZZ=ZZ(1:ns/2+1);
ZZ(2:end-1) = 2*ZZ(2:end-1);
k=ks*(0:ns/2)/ns;
figure,
grid on;
plot(k,YY/2);
hold on
stem(k,ZZ/2);
hold off;
xlim([0 200])
% обнуляем те элменты к спектральной функции...
% к, где спектр полезного сигнала не сильно отличается от нуля, и
% совершаем обратное преобразование Фурье
YY1 = zeros(1,length(k));
for i = 1:length(k)
if k(i) >= 50 && k(i) <= 70
YY1(i) = YY;
elseif k(i) < 50 && k(i) > 70
YY1(i) = 0;
end
end
% обатное преобразование
y = ifft(YY1,ns);
figure
plot(tt,Z,'k-',tt,y+10)
% строим две АКФ
ZZ1 = Z-Y1;
RR = zeros(1,(N-ns-1));
for i = 1:N-ns-1
RR(i) = Z(i)*Z(i+1);
end
RR = (1/(N-i)).*RR;
RR1 = zeros(1,(N-ns-1));
for i = 1:N-ns-1
RR1(i) = ZZ1(i)*ZZ1(i+1);
end
RR1 = (1/(N-i)).*RR1;
figure
plot (i,RR,i,RR1);
clear; close all
TMAX=1;
N=4096;
dt=TMAX/N;
tn=0:dt:TMAX;
A=1;
% Сигнал
% sig=A*(tn>=0.5)+(-A*(tn >= 1));
% sig=A*(tn>=0.5).*(exp(-(5*(tn-0.5).*(tn>=0.5))));
% Сигнал + помеха
x=sig+1*randn(N+1,1)';
L=length(sig); %длина паттерна
N=length(x); %длина входного сигнала
y=zeros(1,N);
for i=1:N
for p=1:L
if (i-p)>0
%y(i)=y(i)+x(i-p)*sig(N-p)*dt;
%y(i)=y(i)+x(p)*sig(N+p-i)*dt;
y(p)=y(i-p)+x(p)*sig(N+p-i)*dt;
end
end
end
% вывод результата
subplot (3,1,1); plot(tn,sig);grid;title('signal') % шаблон сигнала
subplot (3,1,2); plot(tn,x);grid;title('signal+noise') % реальный сигнал
subplot (3,1,3); plot(tn,y);grid;title('filter output') % выход фильтр