2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки


Правила форума


Посмотреть правила форума



Начать новую тему Ответить на тему На страницу 1, 2  След.
 
 Альтернатива оконному Фурье?
Сообщение13.10.2022, 04:40 
Аватара пользователя


26/05/12
1694
приходит весна?
Одна из наиболее часто встречающихся проблем спектрального Фурье-анализа связана с растеканием спектра, когда сигнал не попадает на "гребёнку" базисных функций. Происходит это из-за того, что дискретный Фурье действует на самом деле не на отрезке, а на кольце, получающимся замыканием концов этого отрезка. Когда компонента сигнала не совпадает по частоте с базисной функцией (другими словами, когда длина отрезка разложения не кратна периоду этой компоненты) на стыке концов отрезка у сигнала возникают изломы и разрывы. Как известно, это очень "шумные" особенности в сигнале: разрыв затухает пропорционально $1/f$, а излом (разрыв первой производной) — $1/f^2$. Наглядный пример (зелёный сигнал второго графика — это синий сигнал первого графика, сдвинутый циклически на пол отрезка разложения; на третьем графике — спектры обоих сигналов, наложенные друг на друга):

Изображение

Код в Матлабе, если кто-то захочет поиграться с картинкой:
код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
clc
clearvars
format compact

num = 58;
xx = ((1 : num)' - 0.5) / num;

yy1 = cos (2 * pi * xx * 4.5);
yy2 = [yy1(num / 2 + 1 : end); yy1(1 : num / 2)];

subplot (311)
bar (yy1)
subplot (312)
bar (yy2, 'g')
subplot (313)
bar (abs (fft (yy1)))
hold on
bar (abs (fft (yy2)), 'g')
hold off
 

Решением этой проблемы являются всевозможные оконные функции. Умножение на оконную функцию приводит к тому, что сигнал на стыкуемых концах отрезка зануляется, что устраняет разрывы и изломы. Пики спектров немного уширяются (на две спектральные линии), за то сами пики затухают с удалением от своего центра гораздо быстрее.

Это хорошо проверенное и надёжное решение, которым все пользуются. Однако, у него есть один большой недостаток: потеря информации при занулении участков сигнала. Порой, когда имеющаяся информация и так не велика, такая потеря может быть недопустима. В связи с этим, я подумал: а какая есть альтернатива оконным функциям для борьбы с растеканием спектра? Нельзя ли вместо зануления сигнала, наоборот, дополнить его выдуманным сигналом, чтобы на стыках была непрерывность и гладкость?

На вскидку пришла идея заполнить дополнительную область кольца полиномом, проходящим через стыкуемые концы сигнала. Не очень изящное, но легко реализуемое для эксперимента решение:

Изображение

код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
clc
clearvars
format compact

num = 280;
num_ex = 20;
cyc1 = 5.25;
amp1 = 0.65;
cyc2 = 27.75;
amp2 = -0.3;

nn = (1 : num)';
sgnl =...
    amp1 * sin (2 * pi * cyc1 * nn / num) +...
    amp2 * cos (2 * pi * cyc2 * nn / num);
spct = fft (sgnl);
spct_a = 2 * abs (spct (1 : num / 2 + 1)) / num;
sgnl_ex = [sgnl; zeros(num_ex, 1)];
pnn = [(-4 : 0)' + num, (1 : 5)' + num + num_ex];
pyy = [sgnl((-4 : 0) + num), sgnl(1 : 5)];
pol = polyfit (pnn, pyy, 7);
%sgnl_ex (num + (1 : num_ex)) = sgnl_ex (num) + (sgnl_ex (1) - sgnl_ex (num)) * (1 : num_ex) / num_ex;
sgnl_ex (num + (1 : num_ex)) = polyval (pol, num + (1 : num_ex)');
spct_ex = fft (sgnl_ex);
spct_ex_a = 2 * abs (spct_ex (1 : (num + num_ex) / 2 + 1)) / (num + num_ex);

subplot (2, 1, 1)
plot (sgnl, 'r', 'LineWidth', 2)
hold on
plot ((1 : 100) + num + num_ex, sgnl(1 : 100), 'r', 'LineWidth', 2)
plot ([sgnl_ex; sgnl_ex(1 : 100)], 'k')
hold off
grid on

subplot (2, 1, 2)
plot (0 : num / 2, 20 * log10 (spct_a), 'r', 'LineWidth', 2)
hold on
plot (0 : (num + num_ex) / 2, 20 * log10 (spct_ex_a), 'k')
hold off
ylim ([-100, 0])
grid on
ylabel ('dB')
 


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

Стоит сделать несколько важных наблюдений. Поскольку к сигналу добавляется промежуток с ненулевыми значениями, то суммарная мощность спектра растёт. При этом большая часть спектра теряет мощность. Очевидно, что старая "расплывшаяся" мощность, а так же новая добавленная идут в область с пиками компонент сигнала. Другими словами, спектр стекается обратно в пики. Так же важно отметить, что никакой новой информации в процессе дополнения не возникает. То, как именно дополняется новая область с сигналом, определяется именно самим сигналом, в первую очередь тем, где находятся его пики, то есть, куда мы хотим, чтобы собирался спектр, а где мы хотим, чтобы он занулялся. Это, в некотором смысле, волевое решение, принимаемое на основе имеющегося спектра, таким образом, дополнение сигнала — операция нелинейная. В отличие от оконных функций.

Конечно, используемые мной выше полиномы — это очень неудобные и ненадёжные функции. По многим причинам. Хотелось бы что-то более универсальное, связанное собственно со спектром сигнала. Какие в этом плане могут быть варианты? Что вообще, про такой подход можно сказать? Может, я изобретаю велосипед и всё давно продумано и доказано?

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 09:02 


31/08/22
183
B@R5uk в сообщении #1566609 писал(а):
Нельзя ли вместо зануления сигнала, наоборот, дополнить его выдуманным сигналом, чтобы на стыках была непрерывность и гладкость?

Тоже думал над этим.
Ну или наоборот укоротить выборку до кондиции.
Мысли такие:
Очевидно, что в данной задаче необходимо найти период сигнала. Зная его возможно и вычислить длину на которую нужно дополнить/укоротить выборку до целого периода.
Кроме как пошаговый перебор значения длины придумать не смог.
Можно пробовать автокорреляцию для получения длины.
Понятно, что далеко не все сигналы имеют чисто периодическую структуру и как себя будут вести при этом не периодические сигналы вопрос открытый.

ПС: Если дополнить комплексно сопряженным сигналом той же длины (с разницей на ту же не хватающую длину) что и сама выборка то они сомкнутся в чистое кольцо...
А его можно получить и без вычисления, просто инверсией основной выборки.

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 09:56 


11/08/18
363
B@R5uk в сообщении #1566609 писал(а):
Одна из наиболее часто встречающихся проблем спектрального Фурье-анализа связана с растеканием спектра, когда сигнал не попадает на "гребёнку" базисных функций.

В ЯМРе очень часто возникающая проблема. Обычно решают так: исходный спектр дополняют нулями на $1,...,N$ точек, делают Фурье, скалируют все результаты друг к другу и для каждого пика вытаскивают тот пик, который получается наиболее узким. ИМХО, нулями все-таки не правильно, надо что-то разумное туда добавлять, я собственно тут по соседству как раз это и обсуждаю.

Еще видел "дикий" вариант, когда после Фурье на каждый максимум определяют комплексную экспоненту и фитят эту экспоненту на исходный сигнал в обычных наименьших квадратах. Этот вариант достаточно не стабилен, но, на удивление, большинство пиков удается хорошо зарезолвить.

Кстати, а если сделть вот так:

1. строим Фурье, получаем частоту для какого-то пика. Считаем на сколько у нас период выборки не совпал с целым числом периодов выбранной частоты. Дополняем выборку таким числом отсчетов, но забиваем туда результаты обратного Фурье.
2. повторяем это для каждой частоты, и склеиваем спектры, то есть у нас для каждого пика было свое Фурье, их объединяем и восстанавливаем все назад.

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 10:34 
Заслуженный участник
Аватара пользователя


11/03/08
9905
Москва
Оконные функции с перекрытием. Скажем, Хемминг и перекрытие 50%. Точка, которая войдёт с максимальным весом для другого отрезка войдёт с минимальным. Ничего не теряем (кроме узости пиков).
Дополнять полиномами плохо в принципе, поскольку вносится информация, отсутствующая в исходных данных, и можно произвольно менять ответ. И плохо в плане реализации, поскольку в имеющемся отрезке свободного места для вставки нет, то есть надо произвольно добавить точек, а "гребёнка частот" в результате будет отличаться от считаемой по исходному отрезку.

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 11:28 
Аватара пользователя


26/05/12
1694
приходит весна?
Евгений Машеров в сообщении #1566618 писал(а):
Оконные функции с перекрытием.
Они используются, когда данных наоборот много, и их приходится разбивать на части для анализа. Ситуация совершенно противоположная обсуждаемой.

Евгений Машеров в сообщении #1566618 писал(а):
поскольку вносится информация, отсутствующая в исходных данных, и можно произвольно менять ответ
Это на самом деле не страшно. Произвольность имеет временную локализацию и находится в добавляемой области. После обработки область можно удалить. Даже когда над спектром производятся операции (вырезание/удаление отрезка с пиком, например) Фурье стремится сохранить временную локализацию всех особенностей сигнала (в обратном преобразовании).

Евгений Машеров в сообщении #1566618 писал(а):
то есть надо произвольно добавить точек
Именно! Тут ещё плюс в том, что улучшается частотное разрешение за счёт удлинения кольца разложения. Очень незначительно, конечно, потому что нет особого смысла более чем удваивать длину этого кольца, но всё же.

Евгений Машеров в сообщении #1566618 писал(а):
а "гребёнка частот" в результате будет отличаться от считаемой по исходному отрезку.
Опять же, не страшно. В реальных сигналах данные никогда в точности не попадают на гребёнку, даже если её по одному семплу дополнять нулями, как вон советуют выше. Более того, многие реальные данные имеют значительные флуктуации фазы и амплитуды (не являются чистой синусоидой), поэтому вне зависимости от ухищрений любой пик будет иметь конечную ширину. В пачке спектральных линий и будет зашифрованы эти флуктуации фазы и амплитуды, так что с каким шагом будут идти базисные функции не важно.

Цитата:

Я так понимаю, что никто ничего подобного раньше не придумывал? Uncharted territory?

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 11:44 


11/08/18
363
B@R5uk в сообщении #1566619 писал(а):
Я так понимаю, что никто ничего подобного раньше не придумывал? Uncharted territory?

Так я же написал выше как народ в ЯМРе делает и дополняли с умом спектры, чтобы скачка на было там уже в конце 90-х.

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 11:52 
Аватара пользователя


26/05/12
1694
приходит весна?
ilghiz, а какое название у метода, чтобы загуглить? Есть какие-нибудь образовательные/обзорные книжки/статьи по этому подходу?

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 12:23 


11/08/18
363
B@R5uk в сообщении #1566621 писал(а):
ilghiz, а какое название у метода, чтобы загуглить? Есть какие-нибудь образовательные/обзорные книжки/статьи по этому подходу?

Ох и давно это было... Сам быстро не нагуглил, но очень хорошо помню, что это либо Мартин Биллетер, либо Славка Орехов в Гетеборге мне рассказывали, маловероятно, что они сами эти методы вывели, поэтому и мало вероятно, что в их публикациях это есть, но они возможно пользовали и может быть где-то на эти методы ссылались. Тот же Прони я тоже в той же тусовке услышал, но в ЯМРе его никогда как Прони не называли. Славку постараюсь на днях спросить про этот метод.

Помню, что основная критика тогда этого метода была в большой вычислительной сложности, так как надо было много раз Фурье делать - в 90е годы компьютеры еще не сильно быстрые были. То что я выше предложил

Цитата:
1. строим Фурье, получаем частоту для какого-то пика. Считаем на сколько у нас период выборки не совпал с целым числом периодов выбранной частоты. Дополняем выборку таким числом отсчетов, но забиваем туда результаты обратного Фурье.
2. повторяем это для каждой частоты, и склеиваем спектры, то есть у нас для каждого пика было свое Фурье, их объединяем и восстанавливаем все назад.


по идее должно быть существенно менее трудоемко в случае с не большим числом пиков в спектре - на каждый пик надобно будет только два дополнительных Фурье, одно обратное и оно прямое, использовать.

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 14:44 
Аватара пользователя


26/05/12
1694
приходит весна?
ОК, в результате мозгового штурма получились некоторые продвижения.

Развитие моей идеи заключается в следующем. Вот мы дополнили отрезок с исследуемым сигналом вспомогательным отрезом с подбираемым для непрерывности и гладкости сигналом. По идее, чем длиннее вспомогательный отрезок, тем лучше получится "почистить" спектр исходного сигнала от боковых лепестков, вызванных растеканием. При этом надо будет определиться, какая часть спектра содержит собственно растёкшийся сигнал, а какая должна быть очищена от этого растекания. Варьируя сэмплы вспомогательного отрезка надо подобрать их так, чтобы мощность спектра в области, где сигнал предполагается отсутствующим, была минимальна. Фактически, необходимо зафитить спектром вспомогательного сигнала область спектра с боковыми лепестками основного сигнала и вычесть первое из второго. Всего то!

Результат получается следующим:

Изображение

код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
%   File: study_lacing_2.m
%   Made for: dxdy.ru/post1566622.html
%   Made by: B@R5uk
%   Date: 2022.10.13
clc
clearvars
format compact

num = 280;
num_ex = 70;
num_loop = 50;
cyc1 = 5.25;
amp1 = 0.55;
cyc2 = 27.75;
amp2 = -0.3;
amp_noise = 0.000;

fctr = 1 + num_ex / num;
ppos1 = round (fctr * cyc1);
ppos2 = round (fctr * cyc2);
hwidth = 6;
swidth = floor (num / 2);
swidth_ex = floor ((num + num_ex) / 2);

addr = ones (num + num_ex, 1);
addr ((-hwidth : hwidth) + ppos1 + 1) = 0;
addr ((-hwidth : hwidth) - ppos1 + num + num_ex) = 0;
addr ((-hwidth : hwidth) + ppos2 + 1) = 0;
addr ((-hwidth : hwidth) - ppos2 + num + num_ex) = 0;
addr = 0 ~= addr;

nn = (1 : num)';
sgnl =...
    amp1 * sin (2 * pi * cyc1 * nn / num) +...
    amp2 * cos (2 * pi * cyc2 * nn / num) + ...
    amp_noise * randn (num, 1);
spct = fft (sgnl);
spct_a = 2 * abs (spct (1 : swidth + 1)) / num;

sgnl_ex = [zeros(num_ex, 1); sgnl];
spct_ex = fft (sgnl_ex);
base = [eye(num_ex); zeros(num, num_ex)];
base_spect = fft (base);
fit_vect = base_spect (addr, :) \ spct_ex (addr);
%fit_values = real (ifft (base_spect * fit_vect));
%sgnl_ex = sgnl_ex - fit_values;
%spct_ex = fft (sgnl_ex);
spct_ex = spct_ex - base_spect * fit_vect;
sgnl_ex = real (ifft (spct_ex));
spct_ex_a = 2 * abs (spct_ex (1 : swidth_ex + 1)) / (num + num_ex);

subplot (2, 1, 1)
plot (1 : num_loop, sgnl((1 - num_loop : 0) + num), 'r', 'LineWidth', 2)
hold on
plot (num_loop + (1 : num) + num_ex, sgnl, 'r', 'LineWidth', 2)
plot ([sgnl_ex((1 - num_loop : 0) + num + num_ex); sgnl_ex], 'k')
hold off
grid on

subplot (2, 1, 2)
%plot (20 * log10 (spct_a), 'r', 'LineWidth', 2)
plot ((0 : swidth) / swidth, 20 * log10 (spct_a), 'r', 'LineWidth', 2)
hold on
%plot (20 * log10 (spct_ex_a), 'k')
plot ((0 : swidth_ex) / swidth_ex, 20 * log10 (spct_ex_a), 'k')
plot ((0 : swidth_ex) / swidth_ex, -80 - 20 * addr (1 : swidth_ex + 1), 'b', 'LineWidth', 2)
hold off
ylim ([-100, 0])
grid on
ylabel ('dB')
legend ('Original signal', 'Extended signal', 'Discriminative function')
 


Подавление боковых лепестков на 2 порядка — это, на мой взгляд, успех! И цена не велика: небольшое удлинение отрезка разложения, 2 расчёта БПФ (матрица base_spect содержит комплексные экспоненты, её не обязательно с помощью fft рассчитывать) и обращение матрицы при расчёте фита (самая трудоёмкая часть, в зависимости от длин отрезков). Добавление шума к сигналу процесс не разу не портит, вне зависимости от его величины (устойчивость!). При этом, разумеется, нельзя подавить боковые лепестки ниже уровня шума.

Надо заметить, операция дополнения сигнала всё ещё нелинейная, несмотря на то что все вычисления делаются средствами линейной алгебры. Нелинейность находится в дискриминирующей функции, которая принимает решение о том, какая часть спектра содержит полезный сигнал, а какая — боковые лепестки (в коде эти величины у меня проставляются вручную).

ilghiz в сообщении #1566622 писал(а):
основная критика тогда этого метода была в большой вычислительной сложности, так как надо было много раз Фурье делать
Ещё один плюс только что описанного подхода: от этого недостатка он свободен.

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 17:29 
Заслуженный участник
Аватара пользователя


11/03/08
9905
Москва
Берём отрезок и проводим между начальной и конечной точкой прямую. Вычитая её значения из отсчётов сигнала. После чего начальная и конечная точки равны нулю обе. И разрыва нет.

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 18:21 
Аватара пользователя


26/05/12
1694
приходит весна?
Евгений Машеров в сообщении #1566644 писал(а):
И разрыва нет.
Это у функции нет разрыва. А у её производной? А у второй? Чем выше порядок гладкости, тем быстрей (при удалении от главного максимума) затухают боковые лепестки. Выше я в каком-то смысле (то есть, с фиксированной точностью) достиг совершенства: спектр сигнала локализован вблизи максимумов, соответствующих частотным компонентам, боковых лепестков нет от слова совсем (в пределах заданной точности, конечно), полученная функция, не строго говоря, бесконечно дифференцируема (на сколько это можно утверждать для дискретной функции).

Плюс, если я правильно понял, вы предлагаете добавить сигнал во временной области, соответствующей оригиналу. Это как раз то самое искажение исходного сигнала вспомогательным, когда хотелось бы исследовать оригинал том виде, в каком он имеется изначально. В моём методе добавляемый сигнал локализован вне исходного по построению.

Надо бы ещё проверить на сколько хорошо синусы аппроксимируются в каждом отдельном частотном блоке. Ясно, что при наличии шума, это будет бессмысленный вопрос, но всё же интересна величина интергармонического перекрытия после удаления боковых лепестков.

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 19:42 


11/08/18
363
B@R5uk в сообщении #1566632 писал(а):
Ещё один плюс только что описанного подхода: от этого недостатка он свободен.

И будет у вас 100 пиков и как вы это фиттить будете? Всяко на 100 Фурье попадете.

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 19:58 
Аватара пользователя


26/05/12
1694
приходит весна?
ilghiz, не совсем понял вопрос.

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 20:54 


11/08/18
363
B@R5uk в сообщении #1566621 писал(а):
ilghiz, а какое название у метода, чтобы загуглить? Есть какие-нибудь образовательные/обзорные книжки/статьи по этому подходу?

Начинал делать это Вютрих, это который нобелевку по химии за ЯМР получил, называется peak picking. Последователей много, топталось тоже много, как и публикаций в конце 90-х и чуть позже, но сейчас уже не сильно популярная процедура.

-- 13.10.2022, 19:56 --

B@R5uk в сообщении #1566652 писал(а):
ilghiz, не совсем понял вопрос.

ну как, вы же затачиваете этот метод на одну гармонику (по крайней мере из вашего сообщения я не понял как это сделать на несколько гармоник). Следовательно, если у вас в спектре будет куча пиков, вы всяко на каждый это должны заточить. Так?

 Профиль  
                  
 
 Re: Альтернатива оконному Фурье?
Сообщение13.10.2022, 21:20 
Аватара пользователя


26/05/12
1694
приходит весна?
ilghiz в сообщении #1566657 писал(а):
вы же затачиваете этот метод на одну гармонику
Нет, совсем не так. Сколько угодно гармоник может быть (в разумных пределах), и всё обрабатывается одним махом. На демонстрационных картинках выше две же гармоники. Неужели не заметили?

-- 13.10.2022, 22:04 --

B@R5uk в сообщении #1566647 писал(а):
но всё же интересна величина интергармонического перекрытия после удаления боковых лепестков
Конечно же частотное перекрытие присутствует, чудес не бывает! Поэкспериментировав можно заметить, что чем меньше ширина зануляемого зазора между двумя гармониками, тем сильнее перераспределение мощности между ними (несмотря на более сильное общее подавление лепестков в остальной части спектра). Вот наглядный результат и тестирующая программа:

Изображение

код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
%   File: study_lacing_3.m
%   Made for: dxdy.ru/post1566622.html
%   Made by: B@R5uk
%   Date: 2022.10.13
clc
clearvars
format compact

num = 220;
num_ex = 130;
num_loop = 50;
hwidth = 4;

cyc1 = 5.25;
amp1 = 0.55;
phs1 = 0.0;
cyc2 = 12.75;
amp2 = -0.3;
phs2 = 1.7;
amp_noise = 0.000;

fctr = 1 + num_ex / num;
ppos1 = round (fctr * cyc1);
ppos2 = round (fctr * cyc2);
swidth = floor (num / 2);
swidth_ex = floor ((num + num_ex) / 2);

addr1 = ones (num + num_ex, 1);
addr1 (max (            2,            2 - hwidth + ppos1) : min (swidth_ex - 1,           2 + hwidth + ppos1)) = 0;
addr1 (max (swidth_ex + 1, num + num_ex - hwidth - ppos1) : min (num + num_ex, num + num_ex + hwidth - ppos1)) = 0;
addr1 = 0 ~= addr1;
addr2 = ones (num + num_ex, 1);
addr2 (max (            2,            2 - hwidth + ppos2) : min (swidth_ex - 1,           2 + hwidth + ppos2)) = 0;
addr2 (max (swidth_ex + 1, num + num_ex - hwidth - ppos2) : min (num + num_ex, num + num_ex + hwidth - ppos2)) = 0;
addr2 = 0 ~= addr2;
addr = addr1 & addr2;

addr3 = ones (num, 1);
addr3 (3 : 9) = 0;
addr3 ((-7 : -1) + num) = 0;
addr3 (11 : 19) = 0;
addr3 ((-16 : -9) + num) = 0;

nn = (1 : num)';
sgnl_1 = amp1 * sin (2 * pi * cyc1 * nn / num + phs1);
sgnl_2 = amp2 * sin (2 * pi * cyc2 * nn / num + phs2);
sgnl = sgnl_1 + sgnl_2 + amp_noise * randn (num, 1);
spct = fft (sgnl);
spct_amp = 2 * abs (spct (1 : swidth + 1)) / num;

sgnl_ex = [sgnl; zeros(num_ex, 1)];
spct_ex = fft (sgnl_ex);
base_spect = fft ([zeros(num, num_ex); eye(num_ex)]);
spct_ex_0 = spct_ex - base_spect * (base_spect (addr, :) \ spct_ex (addr));
spct_ex_1 = spct_ex_0;
spct_ex_1 (addr1) = 0;
spct_ex_2 = spct_ex_0;
spct_ex_2 (addr2) = 0;
sgnl_ex_0 = real (ifft (spct_ex_0));
sgnl_ex_1 = real (ifft (spct_ex_1));
sgnl_ex_2 = real (ifft (spct_ex_2));
spct_ex_0_amp = 2 * abs (spct_ex_0 (1 : swidth_ex + 1)) / (num + num_ex);

subplot (6, 1, 1 : 2)
plot ((0 : swidth) / swidth, 20 * log10 (spct_amp), 'g', 'LineWidth', 2)
hold on
plot ((0 : swidth_ex) / swidth_ex, 20 * log10 (spct_ex_0_amp), 'k')
plot ((0 : swidth_ex) / swidth_ex, -100 - 20 * addr (1 : swidth_ex + 1), 'm', 'LineWidth', 2)
hold off
ylim ([-120, 0])
grid on
ylabel ('dB')
legend ('Original signal', 'Extended signal', 'Discriminative function')

subplot (6, 1, 3 : 4)
plot (sgnl, 'g', 'LineWidth', 2)
hold on
plot ([num, num; num + num_ex, num + num_ex]', [-1, 1; -1, 1]', 'k', 'LineWidth', 2)
plot ((1 : num_loop) + num + num_ex, sgnl (1 : num_loop), 'g', 'LineWidth', 2)
plot ([sgnl_ex_1; sgnl_ex_1(1 : num_loop)], 'b')
plot ([sgnl_ex_2; sgnl_ex_2(1 : num_loop)], 'r')
plot ([sgnl_ex_0; sgnl_ex_0(1 : num_loop)], 'k')
hold off
grid on

subplot (6, 1, 5 : 6)
plot (sgnl_ex_1 (1 : num) - sgnl_1, 'b')
hold on
plot (sgnl_ex_2 (1 : num) - sgnl_2, 'r')
plot (0.01 * (real (ifft (spct .* (1 - addr3))) - sgnl), 'g')
hold off
grid on
xlim ([0, num_loop + num + num_ex])
legend ('1st harmonic error', '2nd harmonic error', '1% unextended signal error')
 


Перекрытие проявляется как взаимная противоположность функций ошибок 1-ой и 2-ой гармоник на третьем графике. Сами дополненные гармоники изображены на втором графике синим и красным цветом. Тем не менее, из графиков видно, что ошибка выделения спектра из не дополненного сигнала (зелёная линия на третьем графике) на 2 порядка больше. Конечно, для достижения такого хорошего результата пришлось "нарастить" сигнал более чем на 50%. Так что дополнение свою работу делает: оно действительно собирает "растёкшийся" спектр назад в области локализации пиков компонент сигнала. Даже не смотря на то, что "стягивание" спектра выполняется для всех компонент сразу, то есть для смеси лепестков от всех пиков, каждый пик выбирает и притягивает свою собственную часть растёкшегося спектра. Это, вообще говоря, довольно удивительно.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 25 ]  На страницу 1, 2  След.

Модераторы: Модераторы Математики, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group