Остаточный анализ с автокорреляцией
В этом примере показано, как использовать автокорреляцию с доверительным интервалом, чтобы анализировать остаточные значения, наименьшие квадраты соответствуют к зашумленным данным. Остаточные значения являются различиями между подобранной моделью и данными. В модели шума сигнала-плюс-белый, если у вас есть подходящий вариант для сигнала, остаточные значения должны быть белым шумом.
Создайте набор зашумленных данных, состоящий из полинома 1-го порядка (прямая линия) в аддитивном белом Гауссовом шуме. Аддитивный шум является последовательностью некоррелированых случайных переменных после N (0,1) распределение. Это означает, что все случайные переменные имеют средний нуль и модульное отклонение. Установите генератор случайных чисел на настройки по умолчанию для восстанавливаемых результатов.
x = -3:0.01:3;
rng default
y = 2*x+randn(size(x));
plot(x,y)
Используйте polyfit найти линию наименьших квадратов для зашумленных данных. Отобразите исходные данные на графике наряду с подгонкой наименьших квадратов.
coeffs = polyfit(x,y,1); yfit = coeffs(2)+coeffs(1)*x; plot(x,y) hold on plot(x,yfit,'linewidth',2)
Найдите остаточные значения. Получите последовательность автокорреляции остаточных значений, чтобы отстать 50.
residuals = y - yfit;
[xc,lags] = xcorr(residuals,50,'coeff');Когда вы смотрите последовательность автокорреляции, вы хотите определить, существует ли доказательство автокорреляции. Другими словами, вы хотите определить, похожа ли демонстрационная последовательность автокорреляции на последовательность автокорреляции белого шума. Если последовательность автокорреляции остаточных значений похожа на автокорреляцию белого шумового процесса, вы уверены, что ни один из сигнала не вышел из вашей подгонки и закончился в остаточных значениях. В этом примере используйте 99%-доверительных-интервалов. Чтобы создать доверительный интервал, необходимо знать распределение демонстрационных значений автокорреляции. Также необходимо найти критические значения на соответствующем распределении, между которым находятся 0.99 из вероятности. Поскольку распределение в этом случае является Гауссовым, можно использовать дополнительную обратную функцию ошибок, erfcinv. Отношение между этой функцией и инверсией Гауссовой кумулятивной функции распределения описано на странице с описанием для erfcinv.
Найдите критическое значение для 99%-доверительных-интервалов. Используйте критическое значение, чтобы создать более низкие и верхние доверительные границы.
conf99 = sqrt(2)*erfcinv(2*.01/2); lconf = -conf99/sqrt(length(x)); upconf = conf99/sqrt(length(x));
Постройте последовательность автокорреляции наряду с этими 99%-доверительными-интервалами.
figure stem(lags,xc,'filled') ylim([lconf-0.03 1.05]) hold on plot(lags,lconf*ones(size(lags)),'r','linewidth',2) plot(lags,upconf*ones(size(lags)),'r','linewidth',2) title('Sample Autocorrelation with 99% Confidence Intervals')
Кроме в нулевой задержке, демонстрационные значения автокорреляции лежат в этих 99%-доверительных-границах для автокорреляции белой шумовой последовательности. От этого можно прийти к заключению, что остаточные значения являются белым шумом. А именно, вы не можете отклонить это, остаточные значения являются реализацией белого шумового процесса.
Создайте сигнал, состоящий из синусоиды плюс шум. Данные производятся на уровне 1 кГц. Частота синусоиды составляет 100 Гц. Установите генератор случайных чисел на настройки по умолчанию для восстанавливаемых результатов.
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
rng default
x = cos(2*pi*100*t)+randn(size(t));Используйте дискретное преобразование Фурье (DFT), чтобы получить подгонку наименьших квадратов к синусоиде на уровне 100 Гц. Оценка наименьших квадратов амплитуды равняется 2 / N времена коэффициент ДПФ, соответствующий 100 Гц, где N является длиной сигнала. Действительная часть является амплитудой косинуса на уровне 100 Гц, и мнимая часть является амплитудой синуса на уровне 100 Гц. Подгонка наименьших квадратов является суммой косинуса и синуса с правильной амплитудой. В этом примере интервал ДПФ 101 соответствует 100 Гц.
xdft = fft(x); ampest = 2/length(x)*xdft(101); xfit = real(ampest)*cos(2*pi*100*t)+imag(ampest)*sin(2*pi*100*t); figure plot(t,x) hold on plot(t,xfit,'linewidth',2) axis([0 0.30 -4 4]) xlabel('Seconds') ylabel('Amplitude')
Найдите остаточные значения и определите демонстрационную последовательность автокорреляции, чтобы отстать 50.
residuals = x-xfit;
[xc,lags] = xcorr(residuals,50,'coeff');Постройте последовательность автокорреляции с этими 99%-доверительными-интервалами.
figure stem(lags,xc,'filled') ylim([lconf-0.03 1.05]) hold on plot(lags,lconf*ones(size(lags)),'r','linewidth',2) plot(lags,upconf*ones(size(lags)),'r','linewidth',2) title('Sample Autocorrelation with 99% Confidence Intervals')
Снова, вы видите, что кроме в нулевой задержке, демонстрационные значения автокорреляции лежат в этих 99%-доверительных-границах для автокорреляции белой шумовой последовательности. От этого можно прийти к заключению, что остаточные значения являются белым шумом. А именно, вы не можете отклонить это, остаточные значения являются реализацией белого шумового процесса.
Наконец, добавьте другую синусоиду с частотой 200 Гц и амплитудой 3/4. Соответствуйте только синусоиде на уровне 100 Гц и найдите демонстрационную автокорреляцию остаточных значений.
x = x+3/4*sin(2*pi*200*t);
xdft = fft(x);
ampest = 2/length(x)*xdft(101);
xfit = real(ampest)*cos(2*pi*100*t)+imag(ampest)*sin(2*pi*100*t);
residuals = x-xfit;
[xc,lags] = xcorr(residuals,50,'coeff');Постройте демонстрационную автокорреляцию наряду с этими 99%-доверительными-интервалами.
figure stem(lags,xc,'filled') ylim([lconf-0.12 1.05]) hold on plot(lags,lconf*ones(size(lags)),'r','linewidth',2) plot(lags,upconf*ones(size(lags)),'r','linewidth',2) title('Sample Autocorrelation with 99% Confidence Intervals')
В этом случае значения автокорреляции ясно превышают эти 99%-доверительных-границ для белой шумовой автокорреляции во многих задержках. Здесь можно отклонить гипотезу, что остаточные значения являются белой шумовой последовательностью. Значение - то, что модель не составляла весь сигнал, и поэтому остаточные значения состоят из сигнала плюс шум.