Реализация ФНЧ

@АлександрКом · Регистрация: 21.10.2012

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

Здравствуйте. Задача такая: есть массив данных (запись сигнала). В этом сигнале присутствуют 2 несущие, одна низкочастотная, другая высокочастотная. Необходимо оставить только низкочастотную. Нашел способ, которым можно это сделать - Преобразование Фурье входного сигнала -> преобразование Фурье коэффициентов КИХ-фильтра -> быстрая свертка КИХ-фильтра и сигнала -> обратное преобразование Фурье результатов свертки. Расчет коээфициентов КИХ-фильтра нашел следующий:
должно быть известно заранее
частота среза - Fs
частота дискретизации - Fd
порядок фильтра - n

вычисляем некоторые общие параметры(буквы условные)

w = 2*pi*(Fs/Fd)
q = w/pi

сами коэффициенты, для нечетного порядка расчитываются по ф-ле

А(к) = А(-k) = q*(sin(k*w)/(k*w)), где к = n/2 +/- 1, n/2 +/- 2, n/2 +/- ..., n/2 +/- n/2

где n/2, в данном случае есть ЦЕНТРАЛЬНЫЙ отсчет.

Скажите, будет ли работать данный алгоритм? Может нужно что исправить? Или есть ещё более простой способ?

raxper · 30.03.2014, 10:40

...через накопление и усреднение сигнала тоже можно.

@vital792 · 30.03.2014, 10:49

Сообщение от АлександрКом

Преобразование Фурье входного сигнала -> преобразование Фурье коэффициентов КИХ-фильтра -> быстрая свертка КИХ-фильтра и сигнала -> обратное преобразование Фурье результатов свертки

Это все называется алгоритмом быстрой свертки. То что у вас в середине "быстрая свертка КИХ-фильтра и сигнала" заменяется просто умножением. Вот только при небольшом порядке фильтра, непосредственное выполнение свертки все таки выполняется быстрее этой цепочки. Легко даже посчитать, в каких случаях какой способ быстрее. По поводу вашего фильтра: преобразование фурье вашего фильтра - прямоугольник (если n бесконечно), так что его можно не выполнять, а просто обнулить нужный вам участок. Это простейший ких фильтр с прямоугольным окном, не самый лучший разумеется. Это не раз обсуждалось в данной ветке форума. Из за явления гиббса на границах полосы появляются пульсации. С ними борятся различными способами, самый простой из которых - домножение на весовое окно. Так что советую ничего не изобретать, а просто использовать обычный фильтр и обычную свертку. пример

@eSergey76 · 30.03.2014, 15:06

Более простой способ - WinFilter08.zip. Выбираете тип фильтра и генерите С или VHDL код. Все.

@АлександрКом · 30.03.2014, 16:34 **[ТС]**

Сообщение от vital792

Так что советую ничего не изобретать, а просто использовать обычный фильтр и обычную свертку. пример

Посмотрел пример по ссылке, там вы приводите свой код. Хотел спросить: в коде, я так понял, создается ИХ фильтра и происходит свертка её с сигналом. На выходе получаем отфильтрованный сигнал,да? Хотел узнать по поводу некоторых переменных, правильно ли я понял, что Wn - это частота дискретизации, order - порядок фильтра? Только не понял, наверное, самый главный момент - т.к. это ФНЧ, то должна задаваться частота среза, а в коде как-то не "осознал" этого момента.. Не могли бы пояснить?

@vital792 · 30.03.2014, 17:37

Wn - как раз частота среза. Частота дискретизации тут не важна. Точнее она не используется, т.к. тут нормированная частота, то есть частоту среза задаете от 0 до 1, считая за 1 половину частоты дискретизации

@АлександрКом · 30.03.2014, 20:00 **[ТС]**

Попробовал ваш фильтр при демодуляции ОБП сигнала, что- то не так видно задаю, т.к. на выходе фильтра одни нули. (с Wn - разобрался, но теперь не понятно, что такое lengthFrameForEnvelope

)
Ниже мой код для создания ОБП сигнала с верхней боковой полосой, и его демодуляция:

C++
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
#include <iostream>
#include <cmath>
 
using namespace std;
 
 
#define Wn 0.05 // band for lowpass FIR filter
#define order 20
#define lengthFrameForEnvelope 21
void FIR_Filter(float *x, float *y)
{
    int32_t i, j;
    float b[order+1];//[lengthFrameForEnvelope];
    uint32_t idx = 0;
    float sumB = 0;
    const float piWn = M_PI * Wn;
 
    for(i=-order/2; i<=order/2; i++)
    {
        b[idx] = Wn * sinf(piWn * (i+.00001)) / (piWn * (i+.00001)) *           // Wn * sin(pi*i*Wn) / (pi*i*Wn)
                 (0.54f - 0.46f * cosf(2*M_PI*idx/order));                        // hamming
        sumB += b[idx];
        idx ++;
    }
 
    for(i=0; i<=order; i++)
        b[i] /= sumB;
 
    for(i=0; i<lengthFrameForEnvelope; i++) // convolution
    {
        y[i] = .0;
        for(j=0; j<(i <= order ? i : order); j++)
        {
            y[i] += b[j] * x[i-j];
        }
    }
    for(i=lengthFrameForEnvelope-order/2; i<lengthFrameForEnvelope; i++)
        y[i] = .0;
}
 
 
 
int main()
{
   int Fs=1000;
   float t[1001], s_M[1001], s_SSB_U[1001], y[1001], result[1001];
   t[0]=0;
 
   for(int i=0; i<1000; i++) t[i+1]=t[i]+0.001;
   
int Fc=25, f1=5, f2=10;
  
 for(int j=0; j<1001; j++) {s_M[j]=cos(2*M_PI*f1*t[j])+cos(2*M_PI*f2*t[j]);
   s_SSB_U[j]=cos(2*M_PI*(Fc+f1)*t[j])+cos(2*M_PI*(Fc+f2)*t[j]);
   y[j]=s_SSB_U[j]*cos(2*M_PI*Fc*t[j]);};
 
   FIR_Filter(y, result);
   for(int q=0; q<1001; q++) cout<<s_M[q]<<"       "<<result[q]<<endl;
   return 0;
}

@vital792 · 30.03.2014, 20:41

Сообщение от АлександрКом

(с Wn - разобрался, но теперь не понятно, что такое lengthFrameForEnvelope )

Сообщение от vital792

Добавлено через 1 минуту
lengthFrameForEnvelope вэтом коде - длина выборки сигнала

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

C++
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
#include <iostream>
#include <cmath>
 
using namespace std;
 
#define pi 3.14159265358979
#define Wn 0.05 // band for lowpass FIR filter
#define order 20
#define signalLength 1001
 
void FIR_Filter(float *x, float *y)
{
    int i, j;
    float b[order+1];
    unsigned idx = 0;
    float sumB = 0;
    const float piWn = pi * Wn;
 
    for(i=-order/2; i<=order/2; i++)
    {
        if(i == 0)
            b[idx] = 1.0;
        else
            b[idx] = Wn * sinf(piWn * i) / (piWn * i) *     // Wn * sin(pi*i*Wn) / (pi*i*Wn)
                 (0.54f - 0.46f * cosf(2*pi*idx/order));    // hamming
        sumB += b[idx];
        idx ++;
    }
 
    for(i=0; i<=order; i++)
        b[i] /= sumB;
 
    for(i=0; i<signalLength; i++) // convolution
    {
        y[i] = .0;
        for(j=0; j<(i <= order ? i : order); j++)
        {
            y[i] += b[j] * x[i-j];
        }
    }
    for(i=signalLength-order/2; i<signalLength; i++)
        y[i] = .0;
}
 
int main()
{
   int Fs=1000;
   float t[1001], s_M[1001], s_SSB_U[1001], y[1001], result[1001];
   t[0]=0;
 
   for(int i=0; i<1000; i++) t[i+1]=t[i]+0.001;
   int Fc=25, f1=5, f2=10;
  
 for(int j=0; j<1001; j++)
 {
     s_M[j]=cos(2*pi*f1*t[j])+cos(2*pi*f2*t[j]);
     s_SSB_U[j]=cos(2*pi*(Fc+f1)*t[j])+cos(2*pi*(Fc+f2)*t[j]);
     y[j]=s_SSB_U[j]*cos(2*pi*Fc*t[j]);
 }
 
   FIR_Filter(y, result);
   for(int q=0; q<1001; q++)
       cout<<s_M[q]<<"       "<<result[q]<<endl;
   return 0;
}

@АлександрКом · 30.03.2014, 22:44 **[ТС]**

Спасибо, теперь стало понятнее

Остался последний вопрос: если судить по книжке, по которой я делал свой сигнал, там после демодуляции и отфильтровки сигнал получается такой же, как и тот, который я использовал для модуляции (s_M), только по амплитуде в 2 раза меньше. У меня же в коде, при сравнении значений s_M и result, такой зависимости нет.. В чем может быть проблема?

@vital792 · 30.03.2014, 23:18

Сообщение от АлександрКом

В чем может быть проблема?

в крутизне спада характеристики фильтра. Для такого маленького порядка она довольно плавная и вч колебания фильтруются не полностью.
Вот ваш исходный сигнал:

А вот результат после фильтрации:

Если порядок фильтра увеличить до 40, ачх будет круче и результат будет такой:

@АлександрКом · 01.04.2014, 15:23 **[ТС]**

Спасибо, за объяснение, теперь, действительно стало понятно

Добавлено через 2 часа 7 минут
И снова неприятность

я думал, что порядок фильтра задается #define order 20, но при замене 20 на 40 получил практически тот же результат, что и у вас, но до замены порядка фильтра (когда вч колебания фильтруются не полностью)

@vital792 · 01.04.2014, 17:37

да верно, у меня вышла ошибочка. Вместо

C++
1
2
if(i == 0)
            b[idx] = 1.0;

должно быть

C++
1
2
if(i == 0)
            b[idx] = Wn;

компилятора под рукой не было) При x -> 0 sin(x)/x стремится к 1, но на Wn умножается все выражение.

@АлександрКом · 02.04.2014, 20:04 **[ТС]**

Спасибо большое

@A_Santik · 02.05.2015, 15:08

Реализация полосового фильтра Баттерворта 3 порядка:
ФНЧ-прототип:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?H(S)=\frac{\beta ^{3}}{S^{3}+2\beta S^2+2\beta ^{2}S+\beta^3}$
где
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\beta=\alpha \Omega _{p}= \alpha \cdot tg\left ( \frac{F_s}{F_N} \cdot \frac{\pi }{2}\right )$
Коэффициенты Пф находятся из уравнения, использующего матрицу z-преобразования:
$https://www.cyberforum.ru/cgi-bin/latex.cgi? \begin{bmatrix} 1 &1 & 1 & 1 & 0 &0 & 0\\ 3\Psi & 2\Psi & \Psi &0 & 0&0 &0 \\ 3 &1 &-1 & -3 & 0 & \Psi^2 & 3\Psi^2 \\ 6\Psi & 0 & -2\Psi & 0& 0& 0 & \Psi^3 \\ 3& -1& -1& 3 & 0 & -\Psi^2 &3\Psi^2 \\ 3\Psi & -2\Psi &\Psi & 0 & 0 & 0 &0 \\ -1&1 &-1 & 1 & 0 & 0 & 0 \end{bmatrix}\cdot \begin{bmatrix} 1\\ 2\beta \\ 2\beta^2 \\ \beta^3 \\ 2\beta^2 \\ 2\beta \\ 1 \end{bmatrix}=\begin{bmatrix} A_0\\ A_1\\ A_2\\ A_3\\ A_4\\ A_5\\ A_6 \end{bmatrix}$

где
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\Psi(F_0,F_s) =-2\cdot\frac{cos(\frac{F_0}{F_N}\pi )}{cos(\frac{F_s}{F_N}\frac{\pi }{2})}$
Fo- центральная частота фильтра, Fs - частота среза ФНЧ-прототипа, FN - частота Найквиста.

$https://www.cyberforum.ru/cgi-bin/latex.cgi? [tex]H(z)=\frac{\beta^{3}(1-z^{-2})^3}{A_{0}+A_{1}z^{-1}+A_{2}z^{-2}+A_{3}z^{-3} + A_{4}z^{-4}+A_{5}z^{-5}+A_{6}z^{-6} }$

Т.е. кооэффициенты фильтра можно вычислять прямо в микроконтроллере и оперативно менять частоту/полосу фильтра!

Добавлено через 1 час 13 минут
Если нужен обычный ФНЧ, то расчёты существенно упрощаются:
ФНЧ-прототип:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?H(S)=\frac{\beta ^{3}}{S^{3}+2\beta S^2+2\beta ^{2}S+\beta^3}$
Матрица z-преобразования ФНЧ Баттерворта 3 порядка:
$https://www.cyberforum.ru/cgi-bin/latex.cgi? \begin{bmatrix}1 &1 & 1 & 1 \\ -3 &-1 &1 & 3 \\ 3& -1& -1& 3 \\ -1&1 &-1 & 1\end{bmatrix}\cdot\begin{bmatrix}1\\ 2\beta \\ 2\beta^2 \\ \beta^3\end{bmatrix} = \begin{bmatrix}A_0\\ A_1\\ A_2\\ A_3 \end{bmatrix}$
$https://www.cyberforum.ru/cgi-bin/latex.cgi?H(z)=\frac{\beta^{3}(1+z^{-1})^3}{A_{0}+A_{1}z^{-1}+A_{2}z^{-2}+A_{3}z^{-3} }$

raxper · 02.05.2015, 20:46

...Матлабом не пользуюсь.

А по поводу ресурсов, вы просчитывали, что мало для всех МК? Особенно в тиньке или меге, даже если использовать упрощенные способы из алгоритмических трюков Уоррена. Это вам не ARM с кучей памяти и готовой математикой.

если Вам частоту фильтра надо поменять?

а если не надо? В каких задачах надо менять, никогда не приходилось. Есть Герцель, есть БПФ для случаев анализа всего или части участка.

Я конечно, мог бы формулы окончательные для коэффициентов фильтра

уже есть онлайн-генератор готового Си-кода http://www-users.cs.york.ac.uk... /trad.html

@A_Santik · 04.05.2015, 11:22

ПФ 4-го порядка.
Матрица z-преобразования:
$https://www.cyberforum.ru/cgi-bin/latex.cgi? \begin{bmatrix} 1 & 1 &1 & 1 & 1 & 0 & 0 &0 & 0\\ 4\Psi & 3\Psi & 2\Psi &\Psi &0 & 0 & 0 & 0 &0 \\ 4& 2& 0 & -2 & -4 & 0 & \Psi^2 &3\Psi^2 &6\Psi^2 \\ 12\Psi & 3 & -2\Psi &-3\Psi &0 & 0 &0 &\Psi^3 & 4\Psi^3 \\ 12\Psi^2 &0 &-2 & 0 &6 &0 &-2\Psi^2 &0 &\Psi^4+6 \\ 12\Psi & -3 & -2\Psi &3\Psi &0 & 0 &0 &-\Psi^3 & 4\Psi^3 \\ 4& -2& 0 & 2 & -4 & 0 & \Psi^2 &-3\Psi^2 &6\Psi^2 \\ 4\Psi & -3\Psi &2\Psi & - \Psi & 0& 0& 0 & 0&0 \\ 1& -1& 1 & -1& 1 & 0 & 0 & 0 &0 \end{bmatrix}\begin{bmatrix} 1\\ v\beta \\ w\beta^2 \\ v\beta^3 \\ \beta^4 \\ v\beta^3 \\ w\beta^2 \\ v\beta \\ 1 \end{bmatrix}=\begin{bmatrix} A_0\\ A_1\\ A_2\\ A_3\\ A_4\\ A_5\\ A_6\\ A_7\\ A_8 \end{bmatrix}$
Обозначения, как в предыдущем посте, недостающие коэффициенты в статье http://www.dsplib.ru/content/f... terex.html

Добавлено через 17 часов 5 минут
Матрица z-преобразования для фильтра 15 порядка вот так выглядит:

1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1
-15	-13	-11	-9	-7	-5	-3	-1	1	3	5	7	9	11	13	15
105	77	53	33	17	5	-3	-7	-7	-3	5	17	33	53	77	105
-455	-273	-143	-57	-7	15	17	7	-7	-17	-15	7	57	143	273	455
1365	637	221	21	-43	-35	-3	21	21	-3	-35	-43	21	221	637	1365
-3003	-1001	-143	99	77	-1	-39	-21	21	39	1	-77	-99	143	1001	3003
5005	1001	-143	-187	-11	65	25	-35	-35	25	65	-11	-187	-143	1001	5005
-6435	-429	429	99	-99	-45	45	35	-35	-45	45	99	-99	-429	429	6435
6435	-429	-429	99	99	-45	-45	35	35	-45	-45	99	99	-429	-429	6435
-5005	1001	143	-187	11	65	-25	-35	35	25	-65	-11	187	-143	-1001	5005
3003	-1001	143	99	-77	-1	39	-21	-21	39	-1	-77	99	143	-1001	3003
-1365	637	-221	21	43	-35	3	21	-21	-3	35	-43	-21	221	-637	1365
455	-273	143	-57	7	15	-17	7	7	-17	15	7	-57	143	-273	455
-105	77	-53	33	-17	5	3	-7	7	-3	-5	17	-33	53	-77	105
15	-13	11	-9	7	-5	3	-1	-1	3	-5	7	-9	11	-13	15
-1	1	-1	1	-1	1	-1	1	-1	1	-1	1	-1	1	-1	1

Столбцы матрицы - это коэффициенты произведения полиномов вида (z-1)^n(z+1)^m , где n,m=0...N, n+m=N (порядок фильтра)

Добавлено через 4 часа 46 минут
Внимание!
В матрице z-преобразования ПФ 3 порядка опечатка - изменить знак каждого элемета последней строки.
В матрице z-преобразования ПФ 4 порядка опечатка - 2 столбец, строки 4 и 6: Значения $https://www.cyberforum.ru/cgi-bin/latex.cgi?3\Psi ; -3\Psi$ .

@Red_meat · 14.05.2015, 20:10

Добрый человек, а можно к Вам с еще одним вопросом?
Почему когда Wn = 0.05 это ФНЧ, Wn = 0.5 ФВЧ, а при большинстве прочих выражений это полосной фильтр? И почему это работает, а то что написано у "Айфичер Э. - Цифровая обработка сигналов. Практический подход. 2-е издание. 2004" в методе взвешивания нет?

Добавлено через 21 минуту
Если быть точнее, мне нужно регулировать частоту среза и ширину полосы перехода. Может подскажите как это сделать или что я делаю не так, когда сворачиваю сигнал с фильтром?

@A_Santik · 14.05.2015, 22:25

Я если честно - вопроса не понял совсем

В коэффициентах (для фильтра Баттерворта) очень легко переход ФНЧ-ФВЧ получить. Надо создать ФНЧ с частотой (1-Fs) ( в относительных частотах) и потом поменять знаки у нечётных коэффициентов. С полосовыми фильтрами - сложнее. Могу посоветовать ссылку на работы Костантинидеса.
Приведите конкретный пример фильтра, чтобы я понял суть вопроса

@Red_meat · 14.05.2015, 23:22

Я не про Баттерворта. Я про начала темы и код в ней. Там написан ФНЧ с заранее указанной частотой среза, порядком и прочим. Если хотите - посмотрите.
Что пытаюсь сделать я, генерировать частотные фильтры по типу(ФНЧ, ФВЧ, полосной, заграждающий) частоте среза и затуханию дБ в полосе. Книжка "Айфичер Э. - Цифровая обработка сигналов." говорит о методе взвешивания.
Гласит он, смешать две таблицы(прилагаются) линейной сверткой, и будет вам счастье. Порядок фильтра предлагают брать нечетный, а коэффициенты симметричные. Вроде сказка, а не метод, а получается фигня.

C++
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
//Расчет импульсной характеристики фильтра Fs, Fx края полосы затухания, Fd частота дискретизации
double Fc = (Fs + Fx) / (2 * Fd) ;
for (int i=0;i<=N/2;i++)
{
   if (i==0) H_id[i] = Fc;
   else H_id[i] = sinl(2* M_PI *Fc*i )/(M_PI *i);
    
   // весовая функция Блекмена
   W [i] = 0.42 - 0.5 * cosl((2* M_PI *i) /( N-1)) + 0.08 * cosl((4* M_PI *i) /( N-1));
   H_d [i] = H_id[i] * W[i];
}
 
//Нормировка импульсной характеристики
double SUM=0;
for (int i=0; i<=N/2; i++) SUM +=H_d[i];
for (int i=0; i<=N/2; i++) 
    {
        H_d[i]/=SUM; //сумма коэффициентов равна 1 
        cout<<"H_d[ "<<i<<"] = "<< H_d[i] <<endl;
        if(i!=N/2)
        {
            H[i] = H_d[i];
            H[-i+N-1] = H_d[i];
        }
        else
            H[i] = H_d[i];
        
    }

Если интересно, могу весь код привести, но там дальше просто линейна свертка раскрученная до одного цикла.
Вот. А получается именно то, что описано выше.
Меняешь Fс, и по частотам будто узким окном проходят. При одной это ФНЧ, при другой ФВЧ, а при всех остальных полосной. Но где он окажется, этот не задавленный фильтром участок частот, я не знаю.

И вот или я не очень, или лыжи не едут.

Ну и если уж пошла такая пьянка, может кто подскажет, как АЧХ такого фильтра посмотреть в матлабе можно?

@Red_meat · 16.05.2015, 02:34

Задача решилась! Оставлю код, может кому поможет.

C++
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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
#include "stdafx.h"
 
 
#define _USE_MATH_DEFINES 
#include "stdafx.h"
 
#include <string.h>
#include <iostream>
#include <fstream>
using namespace std;
 
double getWindow(int i, int n, int window)
{
    if( window == 1 )
    {// устраняем нулевые значения
        double a = i - (n - 1) / 2.0;
        if( a < 0 ) a = -a;
        return 2.0 / n * (n / 2.0 - a);
    }
    else if( window == 2 )// устраняем нулевые значения
        return 0.5 - 0.5 * cos(M_PI / n * (1.0 + 2.0 * i));
    if( window == 4 )
    {// устраняем нулевые значения
        double a = M_PI / n * (1.0 + 2.0 * i);
        return 0.5 * (1.0 - 0.16 - cos(a) + 0.16 * cos(2.0 * a));
    }
    else
    {
        double a = 2.0 * M_PI * i / (n - 1);
        if( window == 3 )
            return 0.54 - 0.46 * cos( a );
        else if( window == 5 )
            return 0.35875 - 0.48829 * cos(a) + 0.14128 * cos(2.0 * a) - 0.01168 * cos(3.0 * a);
        else if( window == 6 )
            return 0.35819 - 0.4891775 * cos(a) + 0.1365995 * cos(2.0 * a) - 0.0106411 * cos(3.0 * a);
        else if( window == 7 )
            return 0.355768 - 0.487396 * cos(a) + 0.144232 * cos(2.0 * a) - 0.012604 * cos(3.0 * a);
    }
    //if( type == RECTANGULAR )
    return 1.0;
}
void init(int type, int window, short order, int f1, int f2, int sampleRate,double * m_fir)
{   
 
    if( order == 1 )
    {
        m_fir[0] = 1.0;
        return;
    }
    int n2 = order / 2;
    double w = 2.0 * M_PI * (double)f1 / (double)sampleRate;
    double sum = 0;
    /* расчёт симметричной характеристики для ФНЧ
    double c = (order - 1) / 2.0;
    if( (order&1) != 0 )
    {
    m_fir[n2] = w * getWindow( n2, order, BLACKMAN );
    sum += m_fir[n2];
    }
    for( int i = 0; i < n2; i++ )
    {
    double d = (double)i - c;
    m_fir[i] = Math.sin(w * d) / d * getWindow( i, order, window );
    sum += m_fir[i];
    sum += m_fir[i];
    }
    // нормализация
    if( (order&1) != 0 ) m_fir[n2] /= sum;
    for( int i = 0; i < n2; i++ )
    m_fir[i] = m_fir[order - i - 1] = m_fir[i] / sum;
    */
    int d;
    int i;
    for( i = 0; i < order; i++ )
    {
        d = i - n2;
        if (d !=0)
            m_fir[i] = sin(w * d) / d * getWindow(i, order, window);
        else m_fir[i] = w;
        sum += m_fir[i];
    }
    // нормализация
    for( i = 0; i < order; i++ ) { m_fir[i] /= sum; }
    //
    if( type == 1 ) 
        return; //LOWPASS
    else 
        if( type == 2 ) //HIGHPASS
        {
            for( int i = 0; i < order; i++ ) 
                m_fir[i] = -m_fir[i]; 
 
            m_fir[n2] += 1.0;
            return;
        }
        else
        {// если полосовой или режекторный фильтр
            // расчитываем верхнюю частоту
            double* hf = new double[order];
            w = 2.0 * M_PI * (double)f2 / (double)sampleRate;
            sum = 0;
            for( int i = 0; i < order; i++ )
            {
                int d = i - n2;
                if (d!=0)
                    hf[i] = sin(w * d) / d * getWindow(i, order, window);
                else
                    hf[i] = w;
                sum += hf[i];
            }
            // нормализация
            for( int i = 0; i < order; i++ ) 
                hf[i] /= sum;
            // инвертируем и объединяем с ФНЧ
            for( int i = 0; i < order; i++ ) 
                m_fir[i] -= hf[i];
            m_fir[n2] += 1.0;
            if( type == 3 ) return; //BANDSTOP
            else if( type == 4 ) //BANDPASS
            {
                for( int i = 0; i < order; i++ )
                    m_fir[i] = -m_fir[i]; 
                m_fir[n2] += 1.0;
            }
        }
}
 
int main(int nArgC, char* apszArgV[])
{
    //
    /**
    * инициализация параметров
    * @param type - тип фильтра
    * @param window - окно
    * @param order - порядок фильтра
    * @param f1 - частота ФНЧ и ФВЧ фильтра
    * @param f2 - верхняя частота для полосового и режекторного фильтра
    * @param sampleRate - частота дискретизации
    */
    ifstream in("testin33.txt");
    ofstream out("testout33.txt");
    ofstream HH("H.txt");
 
    int type = 1;
    int window = 2;
    int order = 19;
    int f1 = 1000;
    int f2 = 2000;
    int sampleRate = 8000;
    int sizeIn = 40000;
    
    double *m_fin = new double[order];
    double *H = new double[order];
    double *In = new double[sizeIn];
    double *Out = new double[sizeIn];
    if(order % 2 == 0)
        order -= 1;
    init(type,window,order,f1,f2,sampleRate,H);
 
    int x;
    int i, j;
    for(int c = 0; c < sizeIn; c++)
    {
        in >> x;
        In[c] = x;
 
    }
    for (i=0; i<sizeIn; i++)
    {
        Out[i]=0.;
        for (j=0; (j<order-1) && (j<=i); j++)// та самая формула фильтра
            Out[i]+= H[j]*In[i-j];
    }
    for(int c = 0; c < sizeIn; c++)
    {
        out << Out[c] << " ";
    }//запись в файл
 
    for (int i=0; i<order; i++) 
    {
        HH << H[i] << " ";
        cout<<"H[ "<<i<<"] = "<< H[i] <<endl;
    }
    in.close();
    out.close();
    HH.close();
    cin>>x;
}

Новые блоги и статьи Все статьи Все блоги /
Новый ноутбук volvo 07.12.2025 Всем привет. По скидке в "черную пятницу" взял себе новый ноутбук Lenovo ThinkBook 16 G7 на Амазоне: Ryzen 5 7533HS 64 Gb DDR5 1Tb NVMe 16" Full HD Display Win11 Pro	Музыка, написанная Искусственным Интеллектом volvo 04.12.2025 Всем привет. Некоторое время назад меня заинтересовало, что уже умеет ИИ в плане написания музыки для песен, и, собственно, исполнения этих самых песен. Стихов у нас много, уже вышли 4 книги, еще 3. . .	От async/await к виртуальным потокам в Python IndentationError 23.11.2025 Армин Ронахер поставил под сомнение async/ await. Создатель Flask заявляет: цветные функции - провал, виртуальные потоки - решение. Не threading-динозавры, а новое поколение лёгких потоков. Откат?. . .	Поиск "дружественных имён" СОМ портов Argus19 22.11.2025 Поиск "дружественных имён" СОМ портов На странице: https:/ / norseev. ru/ 2018/ 01/ 04/ comportlist_windows/ нашёл схожую тему. Там приведён код на С++, который показывает только имена СОМ портов, типа,. . .	Сколько Государство потратило денег на меня, обеспечивая инсулином. Programma_Boinc 20.11.2025 Сколько Государство потратило денег на меня, обеспечивая инсулином. Вот решила сделать интересный приблизительный подсчет, сколько государство потратило на меня денег на покупку инсулинов. . . .
Ломающие изменения в C#.NStar Alpha Etyuhibosecyu 20.11.2025 Уже можно не только тестировать, но и пользоваться C#. NStar - писать оконные приложения, содержащие надписи, кнопки, текстовые поля и даже изображения, например, моя игра "Три в ряд" написана на этом. . .	Мысли в слух kumehtar 18.11.2025 Кстати, совсем недавно имел разговор на тему медитаций с людьми. И обнаружил, что они вообще не понимают что такое медитация и зачем она нужна. Самые базовые вещи. Для них это - когда просто люди. . .	Создание Single Page Application на фреймах krapotkin 16.11.2025 Статья исключительно для начинающих. Подходы оригинальностью не блещут. В век Веб все очень привыкли к дизайну Single-Page-Application . Быстренько разберем подход "на фреймах". Мы делаем одну. . .	Фото: Daniel Greenwood kumehtar 13.11.2025	Расскажи мне о Мире, бродяга kumehtar 12.11.2025 — Расскажи мне о Мире, бродяга, Ты же видел моря и метели. Как сменялись короны и стяги, Как эпохи стрелою летели. - Этот мир — это крылья и горы, Снег и пламя, любовь и тревоги, И бескрайние. . .

@АлександрКом 1298 / 927 / 449 Регистрация: 21.10.2012 Сообщений: 2,604

	Реализация ФНЧ 29.03.2014, 23:39. Показов 20884. Ответов 20 Метки нет (Все метки) Здравствуйте. Задача такая: есть массив данных (запись сигнала). В этом сигнале присутствуют 2 несущие, одна низкочастотная, другая высокочастотная. Необходимо оставить только низкочастотную. Нашел способ, которым можно это сделать - Преобразование Фурье входного сигнала -> преобразование Фурье коэффициентов КИХ-фильтра -> быстрая свертка КИХ-фильтра и сигнала -> обратное преобразование Фурье результатов свертки. Расчет коээфициентов КИХ-фильтра нашел следующий: должно быть известно заранее частота среза - Fs частота дискретизации - Fd порядок фильтра - n вычисляем некоторые общие параметры(буквы условные) w = 2pi(Fs/Fd) q = w/pi сами коэффициенты, для нечетного порядка расчитываются по ф-ле А(к) = А(-k) = q(sin(kw)/(k*w)), где к = n/2 +/- 1, n/2 +/- 2, n/2 +/- ..., n/2 +/- n/2 где n/2, в данном случае есть ЦЕНТРАЛЬНЫЙ отсчет. Скажите, будет ли работать данный алгоритм? Может нужно что исправить? Или есть ещё более простой способ? 0

raxper 10236 / 6614 / 498 Регистрация: 28.12.2010 Сообщений: 21,154 Записей в блоге: 1
	30.03.2014, 10:40
	...через накопление и усреднение сигнала тоже можно. 0

@eSergey76 0 / 0 / 0 Регистрация: 26.12.2012 Сообщений: 5
	30.03.2014, 15:06
	Более простой способ - WinFilter08.zip. Выбираете тип фильтра и генерите С или VHDL код. Все. 0

@vital792 2014 / 1286 / 61 Регистрация: 05.06.2010 Сообщений: 2,213
	30.03.2014, 17:37
	Wn - как раз частота среза. Частота дискретизации тут не важна. Точнее она не используется, т.к. тут нормированная частота, то есть частоту среза задаете от 0 до 1, считая за 1 половину частоты дискретизации 0

@АлександрКом 1298 / 927 / 449 Регистрация: 21.10.2012 Сообщений: 2,604
	30.03.2014, 22:44 [ТС]
	Спасибо, теперь стало понятнее Остался последний вопрос: если судить по книжке, по которой я делал свой сигнал, там после демодуляции и отфильтровки сигнал получается такой же, как и тот, который я использовал для модуляции (s_M), только по амплитуде в 2 раза меньше. У меня же в коде, при сравнении значений s_M и result, такой зависимости нет.. В чем может быть проблема? 0

@АлександрКом 1298 / 927 / 449 Регистрация: 21.10.2012 Сообщений: 2,604
	01.04.2014, 15:23 [ТС]
	Спасибо, за объяснение, теперь, действительно стало понятно Добавлено через 2 часа 7 минут И снова неприятность я думал, что порядок фильтра задается #define order 20, но при замене 20 на 40 получил практически тот же результат, что и у вас, но до замены порядка фильтра (когда вч колебания фильтруются не полностью) 0

@АлександрКом 1298 / 927 / 449 Регистрация: 21.10.2012 Сообщений: 2,604
	02.04.2014, 20:04 [ТС]
	Спасибо большое 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 626
	02.05.2015, 15:08
	Реализация полосового фильтра Баттерворта 3 порядка: ФНЧ-прототип: $https://www.cyberforum.ru/cgi-bin/latex.cgi?H(S)=\frac{\beta ^{3}}{S^{3}+2\beta S^2+2\beta ^{2}S+\beta^3}$ где $https://www.cyberforum.ru/cgi-bin/latex.cgi?\beta=\alpha \Omega _{p}= \alpha \cdot tg\left ( \frac{F_s}{F_N} \cdot \frac{\pi }{2}\right )$ Коэффициенты Пф находятся из уравнения, использующего матрицу z-преобразования: $https://www.cyberforum.ru/cgi-bin/latex.cgi? \begin{bmatrix}<br /> 1 &1 & 1 & 1 & 0 &0 & 0\\ <br /> 3\Psi & 2\Psi & \Psi &0 & 0&0 &0 \\ <br /> 3 &1 &-1 & -3 & 0 & \Psi^2 & 3\Psi^2 \\ <br /> 6\Psi & 0 & -2\Psi & 0& 0& 0 & \Psi^3 \\ <br /> 3& -1& -1& 3 & 0 & -\Psi^2 &3\Psi^2 \\ <br /> 3\Psi & -2\Psi &\Psi & 0 & 0 & 0 &0 \\ <br /> -1&1 &-1 & 1 & 0 & 0 & 0<br /> \end{bmatrix}\cdot \begin{bmatrix}<br /> 1\\ <br /> 2\beta \\ <br /> 2\beta^2 \\ <br /> \beta^3 \\ <br /> 2\beta^2 \\ <br /> 2\beta \\ <br /> 1<br /> \end{bmatrix}=\begin{bmatrix}<br /> A_0\\ <br /> A_1\\ <br /> A_2\\ <br /> A_3\\ <br /> A_4\\ <br /> A_5\\ <br /> A_6<br /> \end{bmatrix}$ где $https://www.cyberforum.ru/cgi-bin/latex.cgi?\Psi(F_0,F_s) =-2\cdot\frac{cos(\frac{F_0}{F_N}\pi )}{cos(\frac{F_s}{F_N}\frac{\pi }{2})}$ Fo- центральная частота фильтра, Fs - частота среза ФНЧ-прототипа, FN - частота Найквиста. $https://www.cyberforum.ru/cgi-bin/latex.cgi?<br /> [tex]H(z)=\frac{\beta^{3}(1-z^{-2})^3}{A_{0}+A_{1}z^{-1}+A_{2}z^{-2}+A_{3}z^{-3} + A_{4}z^{-4}+A_{5}z^{-5}+A_{6}z^{-6} }$ Т.е. кооэффициенты фильтра можно вычислять прямо в микроконтроллере и оперативно менять частоту/полосу фильтра! Добавлено через 1 час 13 минут Если нужен обычный ФНЧ, то расчёты существенно упрощаются: ФНЧ-прототип: $https://www.cyberforum.ru/cgi-bin/latex.cgi?H(S)=\frac{\beta ^{3}}{S^{3}+2\beta S^2+2\beta ^{2}S+\beta^3}$ Матрица z-преобразования ФНЧ Баттерворта 3 порядка: $https://www.cyberforum.ru/cgi-bin/latex.cgi? \begin{bmatrix}1 &1 & 1 & 1 \\ -3 &-1 &1 & 3 \\ 3& -1& -1& 3 \\ -1&1 &-1 & 1\end{bmatrix}\cdot\begin{bmatrix}1\\ 2\beta \\ 2\beta^2 \\ \beta^3\end{bmatrix} <br /> = \begin{bmatrix}A_0\\ A_1\\ A_2\\ A_3 \end{bmatrix}$ $https://www.cyberforum.ru/cgi-bin/latex.cgi?H(z)=\frac{\beta^{3}(1+z^{-1})^3}{A_{0}+A_{1}z^{-1}+A_{2}z^{-2}+A_{3}z^{-3} }$ 0

@Red_meat 0 / 0 / 0 Регистрация: 14.05.2015 Сообщений: 9
	14.05.2015, 20:10
	Добрый человек, а можно к Вам с еще одним вопросом? Почему когда Wn = 0.05 это ФНЧ, Wn = 0.5 ФВЧ, а при большинстве прочих выражений это полосной фильтр? И почему это работает, а то что написано у "Айфичер Э. - Цифровая обработка сигналов. Практический подход. 2-е издание. 2004" в методе взвешивания нет? Добавлено через 21 минуту Если быть точнее, мне нужно регулировать частоту среза и ширину полосы перехода. Может подскажите как это сделать или что я делаю не так, когда сворачиваю сигнал с фильтром? 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 626
	14.05.2015, 22:25
	Я если честно - вопроса не понял совсем В коэффициентах (для фильтра Баттерворта) очень легко переход ФНЧ-ФВЧ получить. Надо создать ФНЧ с частотой (1-Fs) ( в относительных частотах) и потом поменять знаки у нечётных коэффициентов. С полосовыми фильтрами - сложнее. Могу посоветовать ссылку на работы Костантинидеса. Приведите конкретный пример фильтра, чтобы я понял суть вопроса 1

Реализация ФНЧ

Решение