clear()
clc()
xdel(winsid()) //zamkniecie okien

fs=1024; //czestotliwosc probkowania
N=1024; //liczba probek
dt=1/fs; //odstep probkowania
T=N*dt; //czas obserwacji sygnalu
t=0:dt:T-dt; //wekot czasu
df=fs/N;
f=0:df:(fs-df);

//parametry syganlow
A1=20;
f1=8;
A2=10;
f2=32;
A3=5;
f3=128;



//Parametry filtru dolnoprzepustowego:
order=6; //rząd filtru
fcutoff=60; //częstotliwość graniczna [Hz]
fcutoffn=(2*fcutoff)/fs; //znormalizowana częstotliwość graniczna [π rad/próbkę]

//Sygnał modulujący
x1=A1*sin(2*%pi*f1*t);
x2=A2*sin(2*%pi*f2*t);
x3=A3*sin(2*%pi*f3*t);
x=x1+x2+x3;
X=fft(x);
fig0=scf();
subplot(311);
plot(t,x);
xlabel('Czas');
ylabel('Amplituda');
title('Sygnał modulujący (informacyjny)');
xgrid;

subplot(312);
plot([f(N/2+1:N)-fs,f(1:N/2)],[abs(X(N/2+1:N)), abs(X(1:N/2))]/N); //przeskalowanie osi
xlabel('Częstotliwość');
ylabel('Amplituda');
title('Widmo amplitudowe dwustronne');
xgrid;

subplot(313);
plot(f(1:32),abs(X(1:32))/fs);
xlabel('Częstotliwość');
ylabel('Amplituda');
title('Fragment widma amplitudowego dwustronnego');
xgrid;

//Sygnał modulowany
y=60*cos(2*%pi*f3*4*t);
Y=fft(y);
fig1=scf();
subplot(311);
plot(t,y);
xlabel('Czas');
ylabel('Amplituda');
title('Sygnał modulowany (nośny)');
xgrid;

subplot(312);
plot([f(N/2+1:N)-fs,f(1:N/2)],[abs(Y(N/2+1:N)), abs(Y(1:N/2))]/N);
xlabel('Częstotliwość');
ylabel('Amplituda');
title('Widmo amplitudowe dwustronne');
xgrid;

subplot(313);
plot(f(48:80),abs(Y(48:80))/fs);
xlabel('Częstotliwość');
ylabel('Amplituda');
title('Fragment widma amplitudowego dwustronnego');
xgrid;

//DWUWSTĘGOWA MODULACJA SYGNAŁU
z=x.*y; // sygnał zmodulowany
Z=fft(z);
fig2=scf();

subplot(311);
plot(t,z);
xlabel('Czas');
ylabel('Amplituda');
title('Sygnał zmodulowany');
xgrid;

subplot(312);
plot([f(N/2+1:N)-fs,f(1:N/2)],[abs(Z(N/2+1:N)), abs(Z(1:N/2))]/N);
xlabel('Częstotliwość');
ylabel('Amplituda');
title('Widmo amplitudowe dwustronne');
xgrid;

subplot(313);
plot(f(40:100),abs(Z(40:100))/fs);
xlabel('Częstotliwość');
ylabel('Amplituda');
title('Fragment widma amplitudowego dwustronnego');
xgrid;

//DEMODULACJA SYGNAŁU
demod=z.*y; //sygnał zdemodulowany

//projektowanie filtru
hz=iir(order,'lp','butt',[fcutoffn/2,0],[.1 .1]); //filtr dolnoprzepustowy Butterwortha
// Odpowiedź częstotliwościowa filtru
[frq,repf]=repfreq(hz,0:0.001:0.5);

fig3=scf();
subplot(211);
plot(t,demod);
xlabel('Czas');
ylabel('Amplituda');
title('Sygnał po demodulacji przed filtracją');
xgrid;

subplot(212);
plot2d(fs*frq,abs(repf),rect=[0,0,fs/2,max(abs(repf))]);
xlabel('Częstotliwość [Hz]');
ylabel('Amplituda');
title('Charakterystyka częstotliwościowa filtru dolnopasmowego');
xgrid;

//filtracja
demodfiltered = flts(demod,hz); //sygnał po demodulacji i filtracji
demodfiltered=2*demodfiltered/(A2^2);
fig4=scf();
subplot(311);
plot(t,demodfiltered);
xlabel('Czas');
ylabel('Amplituda');
title('Sygnał po demodulacji i filtracji');
xgrid;

D=fft(demodfiltered);
subplot(312);
plot([f(N/2+1:N)-fs,f(1:N/2)],[abs(D(N/2+1:N)), abs(D(1:N/2))]/N);
xlabel('Częstotliwość');
ylabel('Amplituda');
title('Widmo amplitudowe dwustronne sygnału po demodulacji');
xgrid;

subplot(313);
plot(f(1:32),abs(D(1:32))/fs);
xlabel('Częstotliwość');
ylabel('Amplituda');
title('Fragment widma amplitudowego dwustronnego');
xgrid;

fig5=scf();
plot(t,x,'r'), plot(t,demodfiltered,'b');
xlabel('Czas');
ylabel('Amplituda');
title('Sygnał modulujący i zdemodulowany');
xgrid;