srnum = 7; % シフトレジスタの数。これにより出力信号の周期が設定される
helz = 2^srnum - 1; % 周期
len = 1000; % 長さ
x = 0:(len-1);

% M系列信号を生成＆表示
msig = maximum_length_sequence(len, srnum);
subplot(3, 1, 1);
plot(x, msig);
title('１．M系列信号');

% M系列信号の自己相関関数を生成＆表示
ac = auto_correlation(0:len, msig);
subplot(3, 1, 2);
plot(x(1:length(ac)), ac);
title('２．自己相関関数');

% M系列の周期毎に縦線を表示（自己相関関数が妥当か、の確認のため）
hold on
    for i = 1 : len

        if helz * i <= (len/2)

            plot([helz*i helz*i] - 1, [min(ac) max(ac)], 'r--'); % x軸の-1は、ズレる分の補正(黒線は'k')
        else

            break;
        end
    end
hold off

% パワースペクトルを表示
subplot(3, 1, 3);
[f, power] = disp_power_spectrum(msig);
plot(f, power);
title('３．パワースペクトル');

function output = maximum_length_sequence(time, srnum)
    % time:時間  srnum:シフトレジスタ数

    % シフトレジスタの数毎の、xorのタップ位置（最後のシフトレジスタと、どれかのシフトレジスタのxor）
    eachtap = [1 0 0 0 0 0 0 0 0;
               1 1 0 0 0 0 0 0 0;
               1 0 1 0 0 0 0 0 0;
               1 0 0 1 0 0 0 0 0;
               0 1 0 0 1 0 0 0 0;
               1 0 0 0 0 1 0 0 0;
               1 0 0 0 0 0 1 0 0;
               0 0 0 0 0 0 0 0 0; % シフトレジスタ数が8の時は、xor1つのみでは（n^2-1周期の信号を生成できない）
               0 0 0 1 0 0 0 0 1];
    z = eachtap(srnum, 1:srnum); % 選択したシフトレジスタの数に合わせたxorのタップ数を抽出
    postap = find(z==1, 1);      % 最後のシフトレジスタ以外のタップ位置を抽出
    output = 0;

    for i = 1 : time

        output(i) = xor(z(postap), z(srnum));

        % それぞれのシフトレジスタの出力を更新
        for j = srnum : -1 : 2

            z(j) = z(j - 1);
        end
        z(1) = output(i); % xor出力と直接繋がるシフトレジスタのみ、xorの出力で初期化
    end

    output(output == 0) = -1; % -1と1の二値信号にするため、0を-1に置換

%     plot(0:length(output) - 1, output); % テスト用

end

function yout = auto_correlation(xin, yin)

    T = length(xin)/2;

    % τを少しずつずらしていく
    for a = 1 : T % a = τ のつもり

        s = 0;

        for t = 2 : T % τを一定値ずらした状態で、2つの波形の積を求める（t = time）

            s = s + yin(t) * yin(t + a); % 右の式を±2/Tの範囲で定積分：x(t)x(t+τ)dt
        end

        yout(a) = s; % x(t)をτの分だけずらした時の積算値を格納
    end

    yout = yout / T;

%     plot(xin(1:length(yout)), yout);
end

function [f, power] = disp_power_spectrum(sig)

    fs = 1; % 帯域
    y = fft(sig);

    n = length(sig);       % サンプル数
    f = (0:n-1) * (fs/n);  % 表示周波数
    power = abs(y) .^ 2/n; % 離散フーリエ変換

    % plot(f,power)
    grid on;
end