Пропустить белый шум через ФНЧ

@Detweeted · Регистрация: 10.10.2012

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

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

Matlab M

F0 = 3e6;                                      % частота среза  
[b, a] = ellip(15, 1, 20, 2*pi*F0, 's');       % синтез коэффициентов фильтра
f = 0:1000:5*10^6;                             % вектор частот для расчета АЧХ и ФЧХ
h = freqs(b, a, 2*pi*f);                       % комплексный коэффициент передачи
subplot(3, 2, 1)
plot(f, abs(h))                                % график АЧХ
grid
subplot(3, 2, 2)
plot(f, unwrap(angle(h))*180/pi);              % график ФЧХ (в градусах)
grid
noise = wgn(1, length(f), -30);                % формируем белый шум
D = fft(h .* fft(noise));                      % пропускаем белый шум через ФНЧ
subplot(3, 2, 3)
plot(f,abs(D));                                % спектр результирующего процесса

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

Зосима · 09.02.2015, 11:25

В строке 12 убери fft

D = h.*fft(noise);

@Detweeted · 18.02.2015, 19:14 **[ТС]**

А вот еще вопрос к знатокам!

А почему если вместо "лобового" обратного преобразования Фурье от спектра (спектр равен h .* fft(noise)) использовать встроенную функцию

D = filter(b, a, noise);

то ничего не будет работать? 0_0

И подскажите еще, пожалуйста. Невероятно мучительный вопрос. А каков математический аналог нормализации корреляционной функции с параметром 'coef'? Ну вот чтоб прям формула была. В хэлпе и интернетах этого нет.

@Detweeted · 22.02.2015, 16:18 **[ТС]**

Ой, Зосимы нет, и никто не отвечает

Сформулирую подробнее.
Почему "лобовой" метод с применением прямого и обратного преобразования фурье работает

Matlab M

F0 = 3e6;                                      % частота среза  
[b, a] = ellip(15, 1, 20, 2*pi*F0, 's');       % синтез коэффициентов фильтра
f = 1:500:5*10^6;                              % вектор частот для расчета АЧХ и ФЧХ
h = freqs(b, a, 2*pi*f);                       % комплексный коэффициент передачи
subplot(3, 2, 1)
plot(f, abs(h))                                % график АЧХ
grid
subplot(3, 2, 2)
plot(f, unwrap(angle(h))*180/pi);              % график ФЧХ (в градусах)
grid
noise = wgn(1, length(f), -30);                % формируем белый шум
D = h .* fft(noise);                           % пропускаем белый шум через ФНЧ
subplot(3, 2, 3)
plot(f,abs(D));                                % спектр результирующего процесса
noise2 = ifft(D);                              % осциллограмма результирующего процесса
subplot(3, 2, 4)
plot(1./f, noise2);                                
xlim([0 10e-3]);

А "универсальный" метод с использованием встроенной функции filter - нет

Matlab M

F0 = 3e6;                                      % частота среза  
[b, a] = ellip(15, 1, 20, 2*pi*F0, 's');       % синтез коэффициентов фильтра
f = 1:500:5*10^6;                              % вектор частот для расчета АЧХ и ФЧХ
h = freqs(b, a, 2*pi*f);                       % комплексный коэффициент передачи
subplot(3, 2, 1)
plot(f, abs(h))                                % график АЧХ
grid
subplot(3, 2, 2)
plot(f, unwrap(angle(h))*180/pi);              % график ФЧХ (в градусах)
grid
noise = wgn(1, length(f), -30);                % формируем белый шум
 % осциллограмма результирующего процесса
noise2 = filter (b, a, noise);
subplot(3, 2, 4)
plot(1./f, noise2);                                
xlim([0 10e-3]);

Зосима · 22.02.2015, 17:14

Detweeted, я твой вопрос видел, но пока не понял в чем тут загвоздка

1. Действительно ли корректна запись:

Matlab M

1
2
3

[b, a] = ellip(15, 1, 20, 2*pi*F0, 's'); 
...
noise2 = filter (b, a, noise);

Может b и a в этих функциях несут разный смысл?

2. А с другим видом ФНЧ (чебышева или батерворта) такая же история? (для данного элиптического в массиве а получается один коэф-т порядка 10¹⁰⁷

)

Зосима · 25.02.2015, 20:42

вот еще один вариант фильтрации:

Matlab M

clear, clc
 
F0 = 3e6;                                      % частота среза  
[b, a] = ellip(15, 1, 20, 2*pi*F0, 's');       % синтез коэффициентов фильтра
f = 1:500:5*10^6;                              % вектор частот для расчета АЧХ и ФЧХ
h = freqs(b, a, 2*pi*f);                       % комплексный коэффициент передачи
subplot(2, 2, 1)
plot(f, abs(h))                                % график АЧХ
grid
subplot(2, 2, 3)
plot(f, unwrap(angle(h))*180/pi);              % график ФЧХ (в градусах)
grid
noise = wgn(1, length(f), -30);                % формируем белый шум
 
% осциллограмма результирующего процесса
H = tf(b, a);
t = linspace(0, 1.5e-4,length(noise));
noise2 = lsim(H,noise,t);
 
subplot(2, 2, 2)
plot(t, noise2,'b')
subplot(2, 2, 4)
plot(t,noise,'r');

@Detweeted · 26.02.2015, 10:59 **[ТС]**

Зосима, выглядит многообещающе, но по факту почему-то не работает

Берем сигнал длительностью 100 микросекунд с несущей 40МГц, пропускаем его через полосовой фильтр с полосой 1МГц. На выходе добрая часть сигнала съедается... Так быть не должно. Не понимаю, откуда такой результат.

Matlab M

clear, clc
Signal = load ('DSignalOVERAL.txt','-ascii');
% частота сигнала  
Fs = 40e6;        
%Длительность сигнала
Ts = 100e-6;
% частота дискретизации
SamplingRate = 10*Fs;
% синтез коэффициентов фильтра
[b, a] = butter(4, [2*pi*(Fs-0.5e6) 2*pi*(Fs+0.5e6)],'s'); 
% вектор частот для расчета АЧХ и ФЧХ
f = 0:1/Ts:SamplingRate;         
% комплексный коэффициент передачи
h = freqs(b, a, 2*pi*f);                       
subplot(2, 2, 1)
plot(f, abs(h))                                % график АЧХ
xlim([Fs-2e6 Fs+2e6]);
grid
subplot(2, 2, 3)
plot(f, unwrap(angle(h))*180/pi);              % график ФЧХ (в градусах)
grid
% осциллограмма результирующего процесса
H = tf(b, a);
t = 0: 1/SamplingRate: Ts;
Signal2 = lsim(H, Signal, t);
 
subplot(2, 2, 2)
plot(t, Signal2,'b')
xlim([0 Ts]);
subplot(2, 2, 4)
plot(t, Signal,'r');
xlim([0 Ts]);

Сам txt-файл с сигналом прикрепил в архиве.

@oleggaaa · 15.04.2015, 22:48

подскажите,пожалуйста, у меня есть исходный сигнал,состоящий из 8 гармоник. Имеется такая схема,состоящая из 4 фильтров Как показать спектральную плотность сигнала,проходящего через каждый фильтр и итоговый сигнал во временной области в конце?

я в матлабе новичок.не очень представляю,как это реализовать.

Зосима 5242 / 3570 / 379 Регистрация: 02.04.2012 Сообщений: 6,473 Записей в блоге: 17
	09.02.2015, 11:25	2
	Сообщение было отмечено Detweeted как решение Решение В строке 12 убери fft *D = h.fft(noise);** 1

@Detweeted 1 / 1 / 0 Регистрация: 10.10.2012 Сообщений: 21
	18.02.2015, 19:14 [ТС]	3
	А вот еще вопрос к знатокам! А почему если вместо "лобового" обратного преобразования Фурье от спектра (спектр равен h .* fft(noise)) использовать встроенную функцию D = filter(b, a, noise); то ничего не будет работать? 0_0 И подскажите еще, пожалуйста. Невероятно мучительный вопрос. А каков математический аналог нормализации корреляционной функции с параметром 'coef'? Ну вот чтоб прям формула была. В хэлпе и интернетах этого нет. 0

@oleggaaa 0 / 0 / 0 Регистрация: 28.02.2015 Сообщений: 30
	15.04.2015, 22:48	8
	подскажите,пожалуйста, у меня есть исходный сигнал,состоящий из 8 гармоник. Имеется такая схема,состоящая из 4 фильтров Как показать спектральную плотность сигнала,проходящего через каждый фильтр и итоговый сигнал во временной области в конце? я в матлабе новичок.не очень представляю,как это реализовать. 0

Пропустить белый шум через ФНЧ

Решение