From ec5005ea97462661ade31a8a4a24143fc3b24b10 Mon Sep 17 00:00:00 2001
From: Unknown <igor.getmanov1@gmail.com>
Date: Mon, 15 Jul 2019 22:37:10 +0300
Subject: [PATCH] =?UTF-8?q?=D0=9D=D0=B0=D0=BF=D0=B8=D1=81=D0=B0=D0=BD=20?=
 =?UTF-8?q?=D0=BA=D0=BE=D0=B4=20=D0=BF=D1=80=D0=BE=D0=B8=D0=B7=D0=B2=D0=BE?=
 =?UTF-8?q?=D0=B4=D1=8F=D1=89=D0=B8=D0=B9=20=D0=BE=D0=BA=D0=BE=D0=BD=D0=BD?=
 =?UTF-8?q?=D0=BE=D0=B5=20=D0=A4=D1=83=D1=80=D1=8C=D0=B5-=D0=BF=D1=80?=
 =?UTF-8?q?=D0=B5=D0=BE=D0=B1=D1=80=D0=B0=D0=B7=D0=BE=D0=B2=D0=B0=D0=BD?=
 =?UTF-8?q?=D0=B8=D0=B5=20(WinFFT.m)=20=D0=B8=20=D0=B2=D0=B8=D0=B7=D1=83?=
 =?UTF-8?q?=D0=B0=D0=BB=D0=B8=D0=B7=D0=B0=D1=86=D0=B8=D1=8E=20(SpectrVisWF?=
 =?UTF-8?q?.m).=20=D0=9C=D0=B5=D1=82=D0=BE=D0=B4=D1=8B=20=D0=B1=D1=8B?=
 =?UTF-8?q?=D0=BB=D0=B8=20=D0=BF=D1=80=D0=BE=D1=82=D0=B5=D1=81=D1=82=D0=B8?=
 =?UTF-8?q?=D1=80=D0=BE=D0=B2=D0=B0=D0=BD=D1=8B=20=D0=B8=20=D0=BF=D1=80?=
 =?UTF-8?q?=D0=BE=D0=B2=D0=B5=D1=80=D0=B5=D0=BD=D1=8B=20=D0=BD=D0=B0=20?=
 =?UTF-8?q?=D0=BE=D1=81=D0=BD=D0=BE=D0=B2=D0=B5=20=D0=B2=D1=81=D1=82=D1=80?=
 =?UTF-8?q?=D0=BE=D0=B5=D0=BD=D0=BD=D1=8B=D1=85=20=D1=84=D1=83=D0=BD=D0=BA?=
 =?UTF-8?q?=D1=86=D0=B8=D0=B9=20Matlab=20(Test1,=20Test2)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 SpectrVisWF.m | 16 +++++++++++++
 Test1.m       | 65 +++++++++++++++++++++++++++++++++++++++++++++++++++
 Test2.m       | 26 +++++++++++++++++++++
 WinFFT.m      | 37 +++++++++++++++++++++++++++++
 4 files changed, 144 insertions(+)
 create mode 100644 SpectrVisWF.m
 create mode 100644 Test1.m
 create mode 100644 Test2.m
 create mode 100644 WinFFT.m

diff --git a/SpectrVisWF.m b/SpectrVisWF.m
new file mode 100644
index 0000000..9a3aa1f
--- /dev/null
+++ b/SpectrVisWF.m
@@ -0,0 +1,16 @@
+function SpectrVisWF(Fsig,t,f)
+%{ 
+    Визуализация спектра в формате "waterfall"
+%}
+
+figure,
+title('Spectrogram')
+ylabel('time [s]')
+xlabel('frequency [Hz]')
+axis([0, Inf,min(t), Inf ])
+hold on
+for ind = 1:size(Fsig,2)
+    plot(f, abs(Fsig(:,ind)) + t(ind),'k')
+    hold on
+end
+hold off
\ No newline at end of file
diff --git a/Test1.m b/Test1.m
new file mode 100644
index 0000000..9f33113
--- /dev/null
+++ b/Test1.m
@@ -0,0 +1,65 @@
+clear, clc, close all
+
+Fs = 8000; % частота квантования
+f1 = 400; % частота заполнения одного из импульсов
+f2 = 1200; % частота заполнения другого импульса
+tau = 0.02; % длительность импульсов
+t1 = 0.02; % время прихода первого импульса
+t2 = 0.06; % время прихода второго импульса
+dt = 1/Fs; % период квантования
+t = 0:dt:0.1; % время
+N = length(t); % количество отсчётов
+% создаём сигналы, состоящие из двух импульсов
+% сигналы отличаются тем, что импульсы с разной частотой
+% приходят в разном порядке
+x1 = (t>=t1 & t<t1+tau).*sin(2*pi*f1*t)+ ...
+ (t>=t2 & t<t2+tau).*sin(2*pi*f2*t);
+x2 = (t>=t1 & t<t1+tau).*sin(2*pi*f2*t)+ ...
+ (t>=t2 & t<t2+tau).*sin(2*pi*f1*t);
+
+% сторим графики сигналов
+figure
+subplot(2,1,1), plot(t,x1,'k'), xlabel('t'), ylabel('x_1(t)')
+subplot(2,1,2), plot(t,x2,'k'), xlabel('t'), ylabel('x_2(t)')
+
+X1 = fft(x1); % дискретное преобразование Фурье первого
+X2 = fft(x2); % дискретное преобразование Фурье второго
+% частоты, соответствующие гармоникам
+f = (ceil(N/2)-N:ceil(N/2)-1)*Fs/N;
+% сторим графики амплитудных спектров сигналов
+figure
+subplot(2,1,1), stem(f, abs(fftshift(X1))/N, 'k.')
+xlabel('f'), ylabel('|X_1|'), axis([0 1600 0 0.12])
+subplot(2,1,2), stem(f, abs(fftshift(X2))/N, 'k.')
+xlabel('f'), ylabel('|X_2|'), axis([0 1600 0 0.12]) 
+% вычисляем оконное дискретное преобразование Фурье
+[Xw1,F,T] = spectrogram(x1,triang(128),64,128,Fs);
+Xw2 = spectrogram(x2,triang(128),64,128,Fs);
+
+c=(1:-0.05:0)'*[1 1 1]; % палитра цветов для построения
+
+% строим контурные графики модуля оконного преобразования
+
+figure
+subplot(2,1,1), contourf(T,F,abs(Xw1))
+colormap(c), colorbar, grid on, axis([min(T) max(T) 0 1600])
+xlabel('t'), ylabel('f'), title('|Xw_1(t,f)|')
+subplot(2,1,2), contourf(T,F,abs(Xw2))
+colormap(c), colorbar, grid on, axis([min(T) max(T) 0 1600])
+xlabel('t'), ylabel('f'), title('|Xw_2(t,f)|')
+
+% вычисляем оконное дискретное преобразование Фурье
+% уменьшили размер окна в 4 раза
+[Xw1,F,T] = spectrogram(x1,triang(32),16,32,Fs);
+Xw2 = spectrogram(x2,triang(32),16,32,Fs);
+
+% строим контурные графики модуля оконного преобразования
+
+% вычисленного с окном другого размера
+figure
+subplot(2,1,1), contourf(T,F,abs(Xw1))
+colormap(c), colorbar, grid on, axis([min(T) max(T) 0 1600])
+xlabel('t'), ylabel('f'), title('|Xw_1(t,f)|')
+subplot(2,1,2), contourf(T,F,abs(Xw2))
+colormap(c), colorbar, grid on, axis([min(T) max(T) 0 1600])
+xlabel('t'), ylabel('f'), title('|Xw_2(t,f)|') 
\ No newline at end of file
diff --git a/Test2.m b/Test2.m
new file mode 100644
index 0000000..cb4ec62
--- /dev/null
+++ b/Test2.m
@@ -0,0 +1,26 @@
+clear all
+
+Fs = 8000; % частота квантования
+f1 = 400; % частота заполнения одного из импульсов
+f2 = 1200; % частота заполнения другого импульса
+tau = 0.02; % длительность импульсов
+t1 = 0.02; % время прихода первого импульса
+t2 = 0.06; % время прихода второго импульса
+dt = 1/Fs; % период квантования
+t = 0:dt:0.1; % время
+N = length(t); % количество отсчётов
+% создаём сигналы, состоящие из двух импульсов
+% сигналы отличаются тем, что импульсы с разной частотой
+% приходят в разном порядке
+x1 = (t>=t1 & t<t1+tau).*sin(2*pi*f1*t)+ ...
+ (t>=t2 & t<t2+tau).*sin(2*pi*f2*t);
+x2 = (t>=t1 & t<t1+tau).*sin(2*pi*f2*t)+ ...
+ (t>=t2 & t<t2+tau).*sin(2*pi*f1*t);
+
+% сторим графики сигналов
+figure
+subplot(2,1,1), plot(t,x1,'k'), xlabel('t'), ylabel('x_1(t)')
+subplot(2,1,2), plot(t,x2,'k'), xlabel('t'), ylabel('x_2(t)') 
+
+[t,f,Fsig] = WinFFT(x1,triang(128), 64, Fs);
+SpectrVisWF(Fsig,t,f)
\ No newline at end of file
diff --git a/WinFFT.m b/WinFFT.m
new file mode 100644
index 0000000..61a72fb
--- /dev/null
+++ b/WinFFT.m
@@ -0,0 +1,37 @@
+function [t,f,Fsig] = WinFFT(sig, wtype, woverlap, Fs)
+%{
+    WindowFFT производит оконное БПФ с параметрами:
+    sig - исходный сигнал;
+    wtype - тип окна, задается в виде массива (рекомендуется выбирать тип
+    окна, равный степени двойки);
+    woverlap - перекрытие окон в отсчетах сигнала;
+    Fs - частота дискретизации сигнала
+
+    Выходные параметры:
+    t - время [c];
+    f - частоты [Гц];
+    Fsig - оконное Фурье-преобразование сигнала sig.
+%}
+
+%{ 
+    Проверка длины сигнала на кратность числу перекрытий. Если длина сигнала
+    не кратна числу перекрытий, то дополняем сигнал нулями справа
+%}
+if mod(length(sig), woverlap) ~= 0                
+    sig = [sig, zeros(1, woverlap * (fix(length(sig)/woverlap) + 1) - length(sig))];
+end
+
+N = length(wtype);                              % количество отсчетов
+f = (0:ceil(N/2) - 1) * Fs /N;                  % частоты, соответствующие гармоникам
+dt = 1/Fs;                                      % период квантования
+k = floor((length(sig) - woverlap)/(N - woverlap));  % число отсчетов по времени с учетом оконного БПФ
+t = (N - woverlap).*dt.*(1:k);                  % время
+
+for ind = 1:k
+    winsig(:,ind) = wtype.*(sig(woverlap * (ind - 1) + 1:woverlap * (ind - 1) + N))';
+end
+
+%% БПФ - берется по столбцам
+Fsig = fft(winsig);
+Fsig = Fsig(1:ceil(N/2),:);
+