2013 / 1285 / 61
Регистрация: 05.06.2010
Сообщений: 2,213
1

краткий обзор методов спектрального оценивания

08.11.2013, 17:33. Показов 5150. Ответов 1
Метки нет (Все метки)

У нас на форуме часто возникают задачи прямо или косвенно связанные со спектральным оцениванием. Решил сделать по этому поводу краткий обзорный топик. Способов оценивания спектра великое множество, всех не описать, поэтому я решил выбрать несколько наиболее распространенных в своем классе.
Что нам нужно от спектрального оценивания, почему этих оценок много и чем они отличаются? Построим например спектр белого шума:
Matlab M
1
2
3
4
x = randn(1, 8196);
spec = fft(x);
psd = abs(spec(1:end/2)).^2;
plot(psd)
получаем нечто жутко изрезанное, не похожее на то что в теории называют спектром белого шума - прямую линию. Поэтому данный способ мало подходит для оценивания спектра подобных сигналов. Если же спектр вычислять окнами с усреднением, то спектр будет сглаживаться и такая спектральная оценка будет более достоверной(метод периодограмм). По сути это то же самое, что пройти по спектру фильтром скользящего среднего, который как известно есть нч фильтр и соответственно сглаживает результат. Однако при этом сглаживаются и острые пики, а значит уменьшается возможность поиска периодического сигнала в шуме.
Первый и самый простой из рассматриваемых методов - непосредственно вычислять преобразование фурье исходного сигнала по конечной сигнальной выборке. Тут думаю пояснения не требуются - этот способ очевиден. Второй практически так же очевиден, когда речь идет о спектральном оценивании случайной последовательности - это преобразование фурье автокорреляционной функции этой последовательности. Если вычислять спектр всего сигнала разом(теоретически бесконечной длины) то первые два метода идентичны. Однако на коротких последовательностях проявляются весьма ощутимые различия, что будет показано ниже. Третий способ относится к классу параметрических методов - авторегрессионный метод (или метод линейного предсказания вперед). В нем подразумевается, что исходный сигнал удовлетворяет некой модели и для оценки спектра требуется оценить параметры этой модели.
Модель представляет из себя белый шум, пропущенный через полюсный фильтр (чисто рекурсивный бих-фильтр). Параметры модели есть коэффициенты этого фильтра (расположены в знаменателе передаточной функции, в числителе только коэффициент усиления)
Оценив коэффициенты этого фильтра, обратным z-преобразованием можно получить нужную нам оценку спектра. Этот метод дает неплохие результаты, если сигнал удовлетворяет указанной модели, например речевой сигнал. И последний из рассмотренных методов - music. Это чисто статистический метод, основанный на разложении корреляционной (в данном случае автокорреляционной) матрицы на собственные числа и вектора. Корреляционная матрица белого шума, как известно диагональна, точнее даже единичная умноженная на дисперсию шума. Если шум будет коррелированным, матрица уже не будет диагональной, однако ее можно привести к диагональному виду разложив на собственные числа и вектора, таким образом отделив коррелированные от некоррелированных, то есть в нашем случае сигнал от шума. По величине собственных чисел (то есть по дисперсиям) можно судить какая часть сигнальная, а какая шумовая и, соответственно выбрать сигнальные и шумовые вектора.
Для моделирования я взял сигнал состоящий из трех синусов, два на близко расположенных частотах и третий подальше, но меньшей амплитуды и белого шума. Ниже код и результаты по выборкам размой длины. В полученных спектрах для простоты опущены всяческие нормировки, то есть результаты отражат качественную, но не количественную оценку спектра. Для наглядности я использовал прямые формулы, без использования быстрых методов (кроме music который и так слишком медленный). Так же отсутствуют весовые окна, без которых на практике не обойтись
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
clear all; close all; clc;
 
k = 2*pi*(0:3999)/8000;
x = sin(k*100)+sin(k*150)+0.5*sin(k*1000); % nF = 51 76 501
% plot(x, 'r'); hold on;
x = awgn(x, -2);
% x = x + 2.0*randn(size(k));
% plot(x);
L = length(x);
L2 = fix(L/2);
p = 100;
rr = xcorr(x,p) ./ p;
 
% periodogram
% nFrames = fix(length(x)/length(rr));
xx = x(1:length(rr));
spm = zeros(1, L2);
for f=0:L2-1
    spm(f+1) = sum(xx(1:length(rr)) .* exp(-2i*pi*(0:length(rr)-1)*f/L));
end
plot(abs(spm.^2))
 
% correlogram
spm = zeros(1, L2);
for f=0:L2-1
    spm(f+1) = sum(rr .* exp(-2i*pi*(0:length(rr)-1)*f/L));
end
figure; plot(abs(spm))
 
% ar
a = levinson(rr(p+1:end),p);
% поиск длины выборки, при которой дисперсия шума в модели перестает расти,
% то есть оптимальной длины. Что то не очень работает)
% var = R(p+1) + sum(a.*R(p+1:end));
% for p=15:5:520
%     R = xcorr(x,p) ./ p;
%     a = levinson(R(p+1:end),p);
%     t = R(p+1) + sum(a.*R(p+1:end));
%     if abs(var - t) < 2
%         break;
%     end
%     var = t;
% end
% var,p
for f=0:L2-1
    tmp = sum(a .* exp(-2i*pi*(0:length(a)-1)*f/L));
    tmp = abs(1+tmp)^2;
    spm(f+1) = 1 / tmp;
end
figure; plot(spm);
 
% music(ev)
nSignals = 3;
R = toeplitz(rr(p+1:end));
[evec, eval] = eig(R);
eval = diag(eval)
noiseEvec = evec(:,nSignals+1:end);
den = 0;
for n = 1:size(noiseEvec, 2)
    h = fft(noiseEvec(:,n), L);
    h = h(1:L2);
    den = den + abs(h).^2./eval(n+nSignals);
end
spm = 1./den;
figure; plot(spm);
Результаты при длине выборки 55. Последовательно расположены спектры полученные методом периодограммы, корелограммы, ар- и music:
краткий обзор методов спектрального оценивания
краткий обзор методов спектрального оценивания
краткий обзор методов спектрального оценивания
краткий обзор методов спектрального оценивания

как видно, только первый метод смог различить две близкие частоты на такой длине выборки, но зато третья гармоника не различима в шуме. У метода music к тому же имеется большая ложная гармоника, не понимаю откуда, видимо в коде где то ошибка, но пока не вижу. Метод хороший, но очень капризный.

Выборка 100:
краткий обзор методов спектрального оценивания
краткий обзор методов спектрального оценивания
краткий обзор методов спектрального оценивания
краткий обзор методов спектрального оценивания

Только второй метод не смог разрешить близкие частоты. В music все та же ложная гармоника, теоретически он должен быть самым лучшим в данных условиях, как будто я где то индексом ошибся и в шумовой базис попал сигнальный вектор

Выборка 200:
краткий обзор методов спектрального оценивания
краткий обзор методов спектрального оценивания
краткий обзор методов спектрального оценивания
краткий обзор методов спектрального оценивания
3
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
08.11.2013, 17:33
Ответы с готовыми решениями:

Краткий обзор Windows Filtering Platform
Здесь я постараюсь дать небольшой обзор WFP, а также затронуть некоторые технические моменты. Без...

Обсуждение "Краткий обзор Windows Filtering Platform"
Статья классная! Она еще актуальна на данный момент? И было бы здорово ссылочку на финальную версию...

Создание системы оценивания
Не могу сделать вывод оценки и количество верных и не верных ответов в тесте.

Блок-схема подпрограммы спектрального анализа
Привет всем! Есть подпрограмма спектрального анализа. function =Spectr_analiz(x, t)...

1
2013 / 1285 / 61
Регистрация: 05.06.2010
Сообщений: 2,213
09.11.2013, 10:20  [ТС] 2
С утра пришел, посмотрел вчерашний код и удивился, как мой music вообще работал? Я же для анализа беру совсем не те вектора. Функция eig() выдает результат сортированный по возрастанию собственных чисел, как и положено для qr алгоритма(правда без использования сдвигов). А в алгоритме используется базис векторов шума, то есть без старших сигнальных. Вот правильная реализация метода music с нормировкой по дисперсиям (ev):
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
nSignals = 3;
R = toeplitz(rr(p+1:end));
[evec, eval] = eig(R);
[eval, idx] = sort(diag(eval), 'descend');
noiseEvec = evec(:, idx);
den = 0;
for n = 1:size(noiseEvec, 2)
    h = fft(noiseEvec(:,n), L);
    h = h(1:L2);
    den = den + abs(h).^2./eval(n);
end
spm = 1./den;
figure; plot(spm);
0
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
09.11.2013, 10:20

Алгоритм рисования спектрального цветового круга
Hello for all. Кто-нибудь знает как можно рисовать спектральный цветовой круг? Меня вполне устроит...

Разработка ПО для спектрального анализа сигналов
Помогите пожалуйста очень нужно разработать программу для разложения в ряд Фурье таблично заданной...

Разработка програми для спектрального анализа
разработка програми для спектрального анализа сигнала, формируемого из совокупности единичних...

написанием функции 100системы оценивания
help! обьясните или покажите хоть как примерно будет выглядить код Написать функцию, которая...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
2
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin
Copyright ©2000 - 2022, CyberForum.ru