Детектирование синусоид в дискретной последовательности

clear all; close all; clc;
 
t=(1:1000)/2000;
f1 = 100; f2 = 200; % частоты
 
s1 = sin(2*pi*f1*t);
s2 = sin(2*pi*f2*t); % блоки сигнала
 
s = [s1 s2 s1 s1 s2 s1 s2]; % сигнал
 
subplot(2,1,1); plot(s); axis tight;
 
fftSize = 1024;         % так как окно не совсем точно совпадает с длиной блока
s = [s zeros(1, 1100)]; % дополняем сигнал нулями
 
% строим спектрограмму
for i=1:floor(length(s)/fftSize)
    tmp = abs(fft(s((i-1)*fftSize+1:i*fftSize)));
    sp(:,i) = tmp(1:end/2); %20*log10(tmp(1:end/2));
end
 
subplot(2,1,2); imagesc(flipud(sp)); grid; axis tight;

clear all; close all; clc;
 
Fs = 2000;
blockLength = 200;
t=(1:blockLength)/Fs;
f1 = 20; f2 = 60; % частоты
 
block(1, :) = sin(2*pi*f1*t);
block(2, :) = sin(2*pi*f2*t); % блоки сигнала
 
% s = [s1 s2 s1 s1 s2 s1 s2]; % сигнал
codeSeq = [0 1 0 0 1 0 1]; % кодовая последовательность, задающая частоты
s = zeros(1, blockLength*length(codeSeq));
stepSignal = blockLength - 50; % перекрытие в 50 отсчетов
for i=1:length(codeSeq)
    s((i-1)*stepSignal+1 : (i-1)*stepSignal+blockLength) = ...
        s((i-1)*stepSignal+1 : (i-1)*stepSignal+blockLength) + block(codeSeq(i)+1, :);
end
s = s ./ max(s);
 
subplot(2,1,1); plot(s); axis tight;
 
fftSize = 128;
 
stepSpectr = 1; % с максимальным перекрытием
% строим спектрограмму
for i=1:floor(length(s)/stepSpectr - fftSize/stepSpectr)
    tmp = abs(fft(s((i-1)*stepSpectr+1:(i-1)*stepSpectr+fftSize)));
    sp(:,i) = tmp(1:end/4); % до Fs/4
end
 
subplot(2,1,2); imagesc(flipud(sp)); grid; axis tight;

function [ ] = Run()
    close all;
    clear all;
    clc;
 
      moon1 = CreateMoon(100, 1, 100, 100);
      moon2 = CreateMoon(100, 27, 64, 50);
      moon3 = CreateMoon(100, 50, 70, 75);
 
      subplot(4, 1, 1);
      plot(moon1, '.');
      axis([0, 100, -2, 2]);
      subplot(4, 1, 2);
      plot(moon2, '.');
      axis([0, 100, -2, 2]);
      subplot(4, 1, 3);
      plot(moon3, '.');
      axis([0, 100, -2, 2]);
      subplot(4, 1, 4);
      plot(moon1+moon2+moon3, '.');
      axis([0, 100, -2, 2]);
  end
 
% function create moon
% dayLength - the number of hours in a day
% appearance - time of appearance
% disappearance - time of disappearance
% period - time for a complete orbit around the planet
function [ y ] = CreateMoon(dayLength, appearance, disappearance, period )
  x = 1:dayLength
  y = zeros(1, length(x));
  y(appearance:disappearance) = y(appearance:disappearance) + sin(x(appearance:disappearance) * 2 * pi / period);
end

moon1 = CreateMoon (1000, 1, 1000, 100);
moon2 = CreateMoon (1000, 340, 1000, 33);
moon3 = CreateMoon (1000, 1, 700, 70);