独自の関数を使用した飽和信号のクリップ除去

ライブスクリプトを開く

データが指定した飽和点より大きい場合、センサーはクリップした読み取り値を返します。読み取り値を再構成するために、飽和区間に隣接する点から多項式をあてはめることができます。再構成を実行する関数を記述し、信号アナライザーに組み込みます。

1 kHz で 14 秒間サンプリングされた 3 チャネルの信号を生成します。信号には、さまざまなサイズと形状の複数のピークがあります。信号を読み取るセンサーは、0.1 V で飽和します。

fs = 1000;
t = 0:1/fs:14-1/fs;

sig = [chirp(t-1,0.1,17,2,'quadratic',1).*sin(2*pi*t/5);
    chirp(t-2,2,2,2.1,'quadratic',100).*exp(-(t-6.5).^2/20).*sin(2*pi*t*2);
    0.85*besselj(0,5*(sin(2*pi*(t+1.5).^2/20).^2)).*sin(2*pi*t/9)]';

sigsat = sig;
stv = 0.1;
sigsat(sigsat >= stv) = stv;

信号アナライザーを開いて、元の信号と飽和した信号を信号テーブルにドラッグします。元のチャネルと飽和したチャネルをそれぞれのディスプレイにドラッグします。

多項式を使用する関数を記述して、信号のピークを再構成します。

最初の入力引数 x は、入力信号です。この引数はベクトルでなければならず、単一チャネルとして取り扱われます。
2 番目の入力引数 tIn は、時間値のベクトルです。ベクトルは信号と同じ長さでなければなりません。入力信号に時間情報がない場合、関数はこの引数を空の配列として読み取ります。
追加の入力引数を指定するには、varargin を使用します。追加の入力引数がない場合は、varargin を省略できます。追加の引数は、[前処理] タブに順序付きのコンマ区切りリストとして入力します。
最初の出力引数 y は、前処理済みの信号です。
2 番目の出力引数 tOut は、出力時間値のベクトルです。入力信号に時間情報がない場合、tOut は空の配列として返されます。
アルゴリズムを実装するために、MATLAB® 関数または Signal Processing Toolbox™ 関数を使用できます。

function [y,tOut] = declip(x,tIn,varargin)
% Declip saturated signal by fitting a polynomial

    % Initialize the output signal

    y = x;

    % For signals with no time information, use sample numbers as abscissas
    
    if isempty(tIn)
        tOut = [];
        t = (1:length(x))';
    else
        t = tIn;
        tOut = t;
    end
    
    % Specify the degree of the polynomial as an optional input argument
    % and provide a default value of 4
    
    if nargin<3
        ndx = 4;
    else
        ndx = varargin{1};
    end

    % To implement your algorithm, you can use any MATLAB or Signal
    % Processing Toolbox function
    
    % Find the intervals where the signal is saturated and generate an 
    % array containing the interval endpoints
    idx = find(x==max(x)); 
    fir = [true;diff(idx)~=1];
    ide = [idx(fir) idx(fir([2:end 1]))];
    % For each interval, fit a polynomial of degree ndx over the ndx+1 points
    % before the interval and the ndx+1 points after the interval
    for k = 1:size(ide,1)
        bef = ide(k,1); aft = ide(k,2);
        intv = [bef-1+(-ndx:0) aft+1+(0:ndx)];
        [pp,~,mu] = polyfit(t(intv),x(intv),ndx);
        y(bef:aft) = polyval(pp,t(bef:aft),[],mu);
    end

end

関数をカスタム前処理関数として信号アナライザーに追加します。[アナライザー] タブの [前処理] ▼ をクリックして、[カスタム関数の追加] を選択します。関数名と説明を入力します。表示されるエディターウィンドウに関数のテキストを貼り付けます。ファイルを保存します。関数が前処理ギャラリーに表示されます。

作成した関数が、飽和した領域を再構成することを確認します。

信号テーブルで、飽和した信号の最初のチャネルを選択します。
[アナライザー] タブの [前処理] ▼ をクリックして、[declip] を選択します。
表示された [前処理] タブで、[前処理] をクリックします。

信号に時間情報がある場合に前処理関数が機能することを確認します。

信号テーブルの sig と sigsat を選択します。個々のチャネルは選択しないでください。
[アナライザー] タブで [時間値] をクリックし、Sample Rate and Start Time を選択して、サンプルレートとして fs を指定します。
上部パネルの信号は、再構成された領域を含み、時間情報を持ちます。