2 / 2 / 0
Регистрация: 08.01.2016
Сообщений: 483
1

Прямое и обратное Фурье преобразование

10.02.2016, 00:32. Показов 5407. Ответов 12
Метки нет (Все метки)

есть исходный сигнал 'w5'. Пусть это звуковой импульс. Известно, что поглощение спектра этого сигнала с расстоянием происходит по следующему закону:
Matlab M
1
H(f,1)=exp(-a.*f.^y.*L(1,n));
где: а-коэффициент поглощения, у-степень (обычно равна 2), f-частота спектра, L- расстояние. То есть поглощение частотно-зависимое.

я пишу код, в котором хочу разложить в спектр исходный сигнал для разных расстояний L с учетом частотно-зависимого поглощения:
Matlab M
1
2
3
4
5
6
7
8
9
10
11
function spec=kerim_sp(w,a,y,L)
% w - вейвлет, 1 столбец
% a - коэф поглощения
% y - степень поглощения
% L - расстояние, одна строка
S=fft(w);
for n=1:length(L);
    f=1:length(S);
    H(f,1)=exp(-a.*f.^y.*L(1,n));
    spec(:,n)=S.*H(f,1);
end
проблема в том, что спектр Фурье отображается симметрично (вроде бы это та же функция, только перевернута относительно оси У) с левой и правой стороны по оси частот.
соответственно если я код в таком виде, в каком он есть, то симметрия нарушается, и восстановить сигнал (обратное преобразование Фурье) не сделать. На рисунке это то, что у меня получается (неправильно получается )

нужно как то изменить код, чтобы функция поглощения H(f,1) применялась к левой и правой стороне спектра симметрично, чтобы затем сделать обратное преобразование Фурье и проследить, как меняется сигнал с расстоянием.

есть предложения, как это реализовать?
Миниатюры
Прямое и обратное Фурье преобразование  
Вложения
Тип файла: rar Фурье.rar (114.6 Кб, 5 просмотров)
__________________
Помощь в написании контрольных, курсовых и дипломных работ, диссертаций здесь
0
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
10.02.2016, 00:32
Ответы с готовыми решениями:

Двумерное прямое и обратное вейвлет-преобразование Хаара, матлаб
Подскажите пожалуйста, как сделать двумерное прямое и обратное вейвлет-преобразование Хаара в...

Обратное преобразование Фурье
Выполняю прямое ПФ, потом прореживаю 0 до размерности обратного ПФ. Проблема в том, что после...

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

Обратное преобразование фурье в матлабе
Цена договорная Требуется взять обратное преобразование фурье от оределенного интеграла, нужен не...

12
Эксперт по электронике
938 / 838 / 121
Регистрация: 23.11.2012
Сообщений: 2,489
10.02.2016, 08:48 2
Лучший ответ Сообщение было отмечено Kerim_Geophysic как решение

Решение

Kerim_Geophysic, дык просто скопируйте одну половину вектора во вторую. Только при этом лучше жестко задать размерность fft таким образом:
Matlab M
1
2
3
4
5
6
7
L = length(w);
nfft = 2^nextpow2(L);
S = fft(w, nfft);
% Потом Ваш код
% ...
% Копирование результатов
spec(:,1:nfft/2) = spec(:,nfft/2:-1:1);
1
529 / 434 / 46
Регистрация: 17.07.2013
Сообщений: 2,223
10.02.2016, 09:02 3
В любом случае строку № 8 из цикла корректнее/надо вынести.
1
2 / 2 / 0
Регистрация: 08.01.2016
Сообщений: 483
10.02.2016, 15:03  [ТС] 4
R2D2, помогла Ваша идея с копированием, я немного переделал, под свою задачу и учел тот факт, что ось симметрии находится на уровне nfft/2+1, в левой части спектра на одну гармонику оказывается больше (вроде нулевая частота, те постоянная составляющая).
вот конечный код:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function spec=kerim_sp(w,a,y,L)
% w - вейвлет, 1 столбец
% a - коэф поглощения
% y - степень поглощения
% L - расстояние, одна строка
x = length(w);
nfft = 2^nextpow2(x);
S = fft(w, nfft);
for n=1:length(L);
    f=1:length(S)
    H(f,1)=exp(-a.*f.^y.*L(1,n));
    spec(:,n)=S.*H(f,1);
end
spec(nfft/2+2:nfft,:) = spec(nfft/2:-1:2,:); % Копирование результатов
end
Nick07, учел Ваш совет, спасибо

возникает другая проблема, не могу сделать обратное преобразование.
тут все очень странно. Если я беру спектр по своему коду для первого столбца ( то есть спектр исходного сигнала, где L=0 и нет затухания) и пытаюсь сделать обратное Фурье преобразование, то у меня получаются комплексные числа.
Matlab M
1
ifft(KomplSpecr(:,1))
Если я буре спектр от исходного сигнала
Matlab M
1
fft(w5)
то я получаю точно такие же числа, в том же порядке, что и в KomplSpecr(:,1), но при этом я могу высчитать обратное Фурье и получить исходный сигнал.
То есть числа одинаковые, но от своего кода я не могу сделать обратное Фурье, а от fft(w5) обратное Фурье получается хорошо. При чем если скопировать данные переменной с fft(w5) в переменную с KomplSpecr(:,1), то обратное Фурье получается.

в чем проблема?
прикрепляю файлы с кодом и рабочим листом
Миниатюры
Прямое и обратное Фурье преобразование  
Вложения
Тип файла: rar Фурье1.rar (85.4 Кб, 5 просмотров)
0
Эксперт по электронике
938 / 838 / 121
Регистрация: 23.11.2012
Сообщений: 2,489
11.02.2016, 01:28 5
Лучший ответ Сообщение было отмечено Kerim_Geophysic как решение

Решение

Цитата Сообщение от Kerim_Geophysic Посмотреть сообщение
то у меня получаются комплексные числа
Лобовое решение:
y = ifft(..., 'symmetric')
А, вообще, конечно, надо делать по уму. И строго проверять симметричность. Тогда и проблем быть не должно.
И кстати, не обязательно гнать цикл. Затухание можно организовать и перемножив Ваш массив спектров и массив в столбцах которого соответственно функции затухания.
1
2 / 2 / 0
Регистрация: 08.01.2016
Сообщений: 483
11.02.2016, 01:50  [ТС] 6
R2D2, спасибо огромное! Вы снова мне очень помогли))
Цитата Сообщение от R2D2 Посмотреть сообщение
И кстати, не обязательно гнать цикл. Затухание можно организовать и перемножив Ваш массив спектров и массив в столбцах которого соответственно функции затухания.
так и сделаю, ато я зациклился на этих циклах.
0
Эксперт по электронике
938 / 838 / 121
Регистрация: 23.11.2012
Сообщений: 2,489
11.02.2016, 02:05 7
Kerim_Geophysic, и правильно)) Циклы в матлабе - зло. Кроме случаев, когда они неизбежны))
fft, кстати, если в нее передать матрицу - вычислит БПФ по столбцам (у Вас это каналы, как я понимаю). Это я так, на всякий случай, вдруг не знали...
1
2 / 2 / 0
Регистрация: 08.01.2016
Сообщений: 483
11.02.2016, 02:12  [ТС] 8
R2D2, и снова спасибо)) потому, что я хотел писать цикл для обратного преобразования по каждому столбцу (как Вы правильно заметили - каналу), но вовремя увидел Ваш ответ)) спасибо за совет))
0
Эксперт по электронике
938 / 838 / 121
Регистрация: 23.11.2012
Сообщений: 2,489
11.02.2016, 17:28 9
Kerim_Geophysic, Ваш файлик с переменными переименовал в spectrums.m (так мне спокойнее). Проверяйте, не исключено, что где то косяк. Просто хотел показать принцип. Понимаю, что код не очевиден, поэтому если с чем то не разберетесь - спрашивайте.
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
clc; close all; clear;
load spectrums
spec = kerim_sp(w5,a,2,L);
a = 1e-5;
n = 2;
nfft = 2^nextpow2(length(w5));
% nfft = length(w5);
Dist = 0:100;
Ls = exp(-a*(1:nfft/2)'.^n*Dist);
Ls = [flipud(Ls); Ls];
S = fft(w5, nfft);
spec1 = (S*ones(1, length(Dist))).*Ls;
%%
figure('Units', 'Normalized', 'Position', [.05 .05 .5 .84])
surf(abs(spec1));
shading interp
colormap summer
camlight(45, 45)
view(-60, 30)
1
Эксперт по электронике
938 / 838 / 121
Регистрация: 23.11.2012
Сообщений: 2,489
11.02.2016, 17:37 10
А на рисунке - результат обратного БПФ. Тоже пришлось использовать параметр 'symmetric'.
Миниатюры
Прямое и обратное Фурье преобразование  
1
2 / 2 / 0
Регистрация: 08.01.2016
Сообщений: 483
11.02.2016, 18:00  [ТС] 11
R2D2,
Цитата Сообщение от R2D2 Посмотреть сообщение
Ваш файлик с переменными переименовал в spectrums.m (так мне спокойнее)
я называю файлы с моим именем, чтобы не перепутать то, что я написал с тем, что я где-то скачал=)) ато я создаю, потом удаляю, если мне не нравится

мне 10 строка непонятна. Я прочитал, что команда flipud(Ls) отобразит симметрично строки переменной Ls. Для чего квадратные скобки там?
0
Эксперт по электронике
938 / 838 / 121
Регистрация: 23.11.2012
Сообщений: 2,489
11.02.2016, 18:31 12
Kerim_Geophysic, квадратные скобки - это так же, как мы задаем обычную матрицу. К примеру
M = [1 2 3 4; 5 6 7 8];
Но здесь не числа в явном виде, а переменные. У нас спектр зеркальный вот и функцию затухания я так же делаю зеркальной. Сначала считаю одну половинку, а потом приставляю к ней вторую но "задом на перед" (если быть точним, то "сверху вниз", т.к. функцию затухания я по столбцам записал).
1
2 / 2 / 0
Регистрация: 08.01.2016
Сообщений: 483
11.02.2016, 18:52  [ТС] 13
R2D2, понял, спасибо
я вчера вот так сделал:
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function spec=kerim_sp_atten(w,a,y,L)
% w - вейвлет, 1 столбец
% a - коэф поглощения
% y - степень поглощения
% L - расстояние, одна строка
x = length(w);
nfft = 2^nextpow2(x);
S = fft(w, nfft);
k=repmat(L,length(S),1);
f=1:length(S)
g=repmat(f',1,length(L));
H=exp(-a.*g.^y.*k);
spec=repmat(S,1,length(L)).*H;
spec(nfft/2+2:nfft,:) = spec(nfft/2:-1:2,:); % Копирование результатов
end
ну, непрофессионально, но сойдет =)
0
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
11.02.2016, 18:52
Помогаю со студенческими работами здесь

Обратное преобразование Фурье, отсекание повторяющегося сигнала
Здравствуйте, Мой вопрос похож на тот, котрый уже тут задавался про затухание сигнала, но на тот...

Как построить графики с нуля по х? (дискретное и обратное преобразования Фурье)
В общем, программа строит мне графики с 1 по оси х. Как сделать, чтобы они строились от 0? Вот...

Дискретное прямое вейвлет преобразование
Здравствуйте, тема такая - Дискретное прямое вейвлет преобразование(дальше мое описание{то как я...

Обратное преобразование
Добрый день! Есть матрица (Создавалась из генератора ПСП) Произвел преобразование по Шаудеру,...


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

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

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