Обратное дискретное преобразование Лапласа

@alexmonik · Регистрация: 05.12.2015

Студворк — интернет-сервис помощи студентам

Здравствуйте, уважаемые пользователи форума.
Долго думал в какой из разделов написать, так как в большей степени вопрос математики, но т.к. реалезуется всё в матлабе, поэтому спрошу здесь.

Сама проблема.
Было получено изображение в области изображений Лапаса вида:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?H(p) = \frac{p}{L{p}^{2}+pR+S}$ ,
L, R, S - действительные коэффициенты,
которое позже будет фигурировать в изображении от двух переменных. Но сначала нужно разобраться с ним.
Из данного изображения необходимо получить оригинал функции(разумеется, в дискретной форме; получать оригинал в аналитическом виде нет необходимости, да и в дальнейшем будет не целесообразно).
Т.к., в MATLAB нет готового обратного дискретного преобразования Лапласа, я так понимаю, необходимо произвести замену
$https://www.cyberforum.ru/cgi-bin/latex.cgi?p = j\omega$ ,
где j - мнимая единица, $https://www.cyberforum.ru/cgi-bin/latex.cgi?\omega$ - частота; и после этого воспользоваться обратным быстрым преобразоваением Фурье ifft(), чтобы получить оригинал изображения.
Разумеется, выполнив в лоб замену и применяя ifft(), получается ерунда.
Тренировался на злощастной синусоиде, но и с ней ничего не вышло. Далее привожу алгоритм, код и результаты.
1. взял изображение sin(x)
$https://www.cyberforum.ru/cgi-bin/latex.cgi?L{[sin(X)]} = \frac{1}{p^2+1}$

2. далее строю АЧХ (для проверки тот ли спектр получил)

Matlab M
1
2
3
4
5
6
7
f = 500; %   частота сигнала
j = sqrt(-1);% мнимая единица
Fd = 2048;  % частота дискретизации
FftL = Fd*4;    % количество линий спектра
omeg = -Fd/2:Fd/FftL:Fd/2-Fd/FftL; % частота
H = f./((j*omeg).^2+f^2); % изображение Лапласа, с заменой p = j*omeg
plot(omeg,abs(H)/2)

Обратное дискретное преобразование Лапласа

вроде всё верно

3. далее для ifft нужно произвести манипуляции со спектром, чтобы нулевая частота находилась в начале вектора частот

Matlab M
1
2
3
h1(1:FftL/2+1) = H(FftL/2:FftL); % манипуляции с отзеркаливанием спектра
h1(FftL/2+1:FftL) = H(1:FftL/2); %  ---//---
s = ifft(H,FftL);   % обратное преобразование

и вот здесь начинаются проблемы, а именно все значения s = NaN, т.е, где-то, как я понимаю, выскакиевает деление на ноль.

Буду благодарен, если ткенёте меня носом в ошибку.

и еще раз весь скрипт целиком:

Matlab M
1
2
3
4
5
6
7
8
9
10
f = 500; %   частота сигнала
j = sqrt(-1);% мнимая единица
Fd = 2048;  % частота дискретизации
FftL = Fd*4;    % количество линий спектра
omeg = -Fd/2:Fd/FftL:Fd/2-Fd/FftL; % частота
H = f./((j*omeg).^2+f^2); % изображение Лапласа, с заменой p = j*omeg
h1(1:FftL/2+1) = H(FftL/2:FftL); % манипуляции с отзеркаливанием спектра
h1(FftL/2+1:FftL) = H(1:FftL/2); %  ---//---
s = ifft(H,FftL);   % обратное преобразование
plot(omeg,abs(H)/2)

@alexmonik · 28.02.2018, 18:38 **[ТС]**

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

Matlab M
1
2
s = ifft(h1,FftL);   % обратное преобразование
plot(s)

не исправил когда спектр строил.

@Annza · 01.03.2018, 13:28

Есть тонкости с этим отзеркаливанием спектра. Они связаны с тем, что ДПФ в Матлабе для четного количества отсчетов, например N = 8, выдаст вектор. 1ая компонента вектора - постоянная составляющая сигнала, со 2ой по N/2 компоненту - часть спектра в положительной области частот, N/2+1 - самосопряженный отсчет(всегда действительный) на частотах +-Fs/2, с N/2+2 по N компоненту - симметричная отрицательная часть спектра. Вот тут пытались разобраться

В вашем случае в h1 получается Inf - потому что в знаменателе L при некоторых значениях 0. Это так действительно по теории? Я просто не знаю

В общем я тут по-другому отзеркаливаю и Inf ограничиваю значением 1000:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
clear all
close all
f = 500; %   частота сигнала
j = sqrt(-1);% мнимая единица
Fd = 2048;  % частота дискретизации
FftL = Fd*4;    % количество линий спектра
omeg = -Fd/2:Fd/FftL:Fd/2-Fd/FftL; % частота
H = f./((j*omeg).^2+f^2); % изображение Лапласа, с заменой p = j*omeg
 
H = f./((j*omeg).^2+f^2); % изображение Лапласа, с заменой p = j*omeg
figure, plot(abs(H)/2)
 
h1(1:FftL/2) = H(FftL/2+1:FftL); % манипуляции с отзеркаливанием спектра
h1(FftL/2+1:FftL) = H(1:FftL/2); %  ---//---
figure, plot(abs(h1)/2)
h1(h1==Inf)=1000;
 
s = ifft(h1,FftL);   % обратное преобразование
figure, plot(s)

@alexmonik · 15.05.2018, 11:05 **[ТС]**

Итак, спустя почти 3 месяца

я решил свою проблему. Поэтому, если у кого-то есть аналогичная проблема, а именно получение из изображения в области Лапласа оригинала с помощью обратного дискретного преобразования Фурье, то вот как я это делал. (это сообщение для таких же глупеньких как я, так что за пояснение совсем элементарных вещей не пинайте)
к примеру, у нас есть сигнал вида $https://www.cyberforum.ru/cgi-bin/latex.cgi?cos(\omega t)$ длительностью [0;T] и у нас есть его изображение в области Лапласа. Теперь, из этого изображения мы хотим получить обратно косинус, но не через обратное преобразование Лапласа, а через обратное дискретное преобразования Фурье.
Первое, что стоит отметить, что табличные изображения даются на случай сигнала бесконечной длительности, поэтому:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\int_{0}^{\infty }cos(ft)exp(-st)dt = \frac{s}{s^2+f^2}$ ,
а в случае конечной длительности уже:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\int_{0}^{T }cos(ft)exp(-st)dt = \frac{s}{s^2+f^2} + \frac{fsin(fT)-scos(fT)}{s^2+f^2}exp(-sT)$
Таким образом мы получили изображение нашего сигнала. Теперь чтобы можно было применить ОДПФ, перейдём от изображения Лапласа к фурье-образу. Для чётных функций справедлива замена: $https://www.cyberforum.ru/cgi-bin/latex.cgi?s = 2\pi j\omega$ .
омега, которая пришла из подстановки $https://www.cyberforum.ru/cgi-bin/latex.cgi?s = 2\pi j\omega$ будет переменной величиной. С ней-то у меня и было больше всего проблем. Строго говоря, конечный сигнал имеет бесконечный спектр, но т.к. у нас нет бесконечного количества времени, частоту нужно брать:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\omega = [0:\frac{1}{T}:kf+1], k\geq 2$
(почему так, можно посмотреть, к примеру, в [Макс Ж. - Методы и техника обработки сигналов при физических измерениях] и [Сергиенко А.Б. - Цифровая обработка сигнала] )
И, собственно всё, можно напрямую подставлять наше изображение в функцию ifft(). Ниже код.

Matlab M
1
2
3
4
5
6
7
8
9
10
11
j1 = sqrt(-1); % мнимая единица
T = 6*pi; % длительность сигнала
f = 1; % частота сигнала
omeg = 0:1/T:20*f+1; % вектор частот
p = 2*pi*j1*omeg; % замена переменной p на j*omega
H = (p./(p.^2+f^2)) + (exp(-p*T).*(f.*sin(f*T)-p.*cos(f*T)))./(p.^2+f^2); % изображение сигнала
H(~isfinite(H))= 1; % может случиться, что числитель в верхней строке будет равен нулю
t = 0:T/(length(omeg)-1):T; % вектор времен для построения оригинала изображения
plot(t,real(ifft(H*length(omeg)))) % во втором аргументе функции plot() сразу применяется функция обратного преобразования 
                                                 %фурье ifft() и массив с изображением домножается на собственную длину (т.к. в ifft() 
                                                 %производится деление на эту самую длину)

Новые блоги и статьи Все статьи Все блоги /
Мысли в слух. Про "навсегда". kumehtar 16.04.2026 Подумалось тут, что наверное очень глупо использовать во всяких своих установках понятие "навсегда". Это очень сильное понятие, и я только начинаю понимать край его смысла, не смотря на то что давно. . .	My Business CRM MaGz GoLd 16.04.2026 Всем привет, недавно возникла потребность создать CRM, для личных нужд. Собственно программа предоставляет из себя базу данных клиентов, в которой можно фиксировать звонки, стадии сделки, а также. . .	Знаешь почему 90% людей редко бывают счастливыми? kumehtar 14.04.2026 Потому что они ждут. Ждут выходных, ждут отпуска, ждут удачного момента. . . а удачный момент так и не приходит.	Фиксация колонок в отчете СКД Maks 14.04.2026 Фиксация колонок в СКД отчета типа Таблица. Задача: зафиксировать три левых колонки в отчете. Процедура ПриКомпоновкеРезультата(ДокументРезультат, ДанныеРасшифровки, СтандартнаяОбработка) / / . . .
Настройки VS Code Loafer 13.04.2026 { "cmake. configureOnOpen": false, "diffEditor. ignoreTrimWhitespace": true, "editor. guides. bracketPairs": "active", "extensions. ignoreRecommendations": true, . . .	Оптимизация кода на разграничение прав доступа к элементам формы Maks 13.04.2026 Алгоритм из решения ниже реализован на нетиповом документе, разработанного в конфигурации КА2. Задачи, как таковой, поставлено не было, проделанное ниже исключительно моя инициатива. Было так:. . .	Контроль заполнения и очистка дат в зависимости от значения перечислений Maks 12.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "ПланированиеПерсонала", разработанного в конфигурации КА2. Задача: реализовать контроль корректности заполнения дат назначения. . .	Архитектура слоя интернета для сервера-слоя. Hrethgir 11.04.2026 В продолжение https:/ / www. cyberforum. ru/ blogs/ 223907/ 10860. html Знаешь что я подумал? Раз мы все источники пишем в голове ветки, то ничего не мешает добавить в голову такой источник, который сам. . .

Обратное дискретное преобразование Лапласа

Решение