Фильтрация

seq = [1 0 1 0 1 0 ]; % некая кодовая последовательность(в данном случае М)
numberOfPeriods = 2; % длина элементарного импульса в периодах несущей
phaseSteps = seq * pi;
 
Ampl = 5; % ам не используется
 
Fs = 500000; % частота дискретизации
fc = 100; % несущая
 
chipDuration = numberOfPeriods / fc; % длительность элементарного импульса
fullLength = ceil(length(seq) * chipDuration * Fs); % длина сигнала в отсчетах
phase = zeros(1, fullLength);
 
 
 
phase = zeros(1, fullLength);
currPhase = phaseSteps(1);
FlipIdxReal = chipDuration * Fs;
 
ChipIdx = 1;
PrevChipIdx = 1;
 
for i=1:fullLength
   
    phase(i) = currPhase;
    
    currPhase = currPhase + 2*pi*fc/Fs;
    if i>FlipIdxReal
        PrevChipIdx = ChipIdx;
        ChipIdx = ChipIdx + 1;
     
        currPhase = currPhase + phaseSteps(ChipIdx) - phaseSteps(PrevChipIdx);
        FlipIdxReal = FlipIdxReal + chipDuration * Fs;
    end;
end;
 
signal = Ampl * sin(phase);
noiseLevel = 1;
noise = noiseLevel * randn(1, length(signal));
noisedSignal = signal + noise;
 
subplot(2,1,1);
 
[a, b]=size(seq);
seq(b+1) = 1;
stairs(seq) ;
AXIS([1 7 -1 2])
title ('Исходный сигнал');
subplot(2,1,2);
plot(noisedSignal);
title ('Фазовая манипуляция');
 
figure;
phase(1) = 1;
for i=2:fullLength
   
    phase(i) = phase(i-1) + 2*pi*fc/Fs;  
end
sample1=Ampl * sin(phase)
d = (sample1 - noisedSignal).^2;
plot(d);
 
k=1;
 
FlipIdxReal=FlipIdxReal-FlipIdxReal;
for i=1:fullLength
    if i>FlipIdxReal
        if d(i)>50
            nseq(k)= 1
        else 
            nseq(k)=0
        end;
        k = k+1;
        FlipIdxReal = FlipIdxReal + chipDuration * Fs;
    end;
end;
 
f = 0 : Fs/fullLength : Fs - Fs/fullLength;
figure
xw=abs(fft(d));
plot(f(1 : fullLength/2), xw(1 : fullLength/2))
title('Спектор')
 
 
 
a0=1;
a1=5.5610;
a2= 0.6414;
 
b0=0.0201;
b1=0.0402;
b2=0.0201;
 
A=[a0, a1, a2];
B=[b0, b1, b2];
f1 = 100;
figure;
y1 = filter(B, A, d);
%y2 = ifft(fft(d).*fft(B, fullLength)./fft(A, fullLength));
plot(y1);
axis([0 2/f1 -1 1]);
title('Filter in time')
 
figure;
spy1 = abs(fft(xw));
plot(f, spy1, 'g');
title('Filter in frequency')
 
inT = ifft(spy1);
figure;
plot (inT);
AXIS([1 70 0 100])

phase(1) = 1;
for i=2:fullLength
   
    phase(i) = phase(i-1) + 2*pi*fc/Fs;  
end
sample1=Ampl * sin(phase)
d = (sample1 - noisedSignal).^2;
plot(d);

FlipIdxReal=FlipIdxReal-FlipIdxReal;
for i=1:fullLength
    if i>FlipIdxReal
        if d(i)>50
            nseq(k)= 1
        else 
            nseq(k)=0
        end;
        k = k+1;
        FlipIdxReal = FlipIdxReal + chipDuration * Fs;
    end;
end;

...
title ('Фазовая манипуляция');
 
ss = noisedSignal .* sin(2*pi*fc*(0:length(noisedSignal)-1)/Fs);
% figure; plot(ss)
FsFc = Fs/fc;
b = ones(1, FsFc)/FsFc;
envelope = filter(b, 1, ss);
% figure; plot(envelope)
chipLen = chipDuration * Fs;
k=1;
for i=0:fullLength-1
    if mod(i, chipLen) == 0
        nseq(k) = mean(envelope(i+1:i+chipLen));
        k = k+1;
    end
end
nseq = (1-sign(nseq))/2;
figure; stairs([nseq 0]); axis([1 length(nseq)+1 -1 2])