Форум программистов, компьютерный форум, киберфорум
Наши страницы
Цифровая обработка сигналов
Войти
Регистрация
Восстановить пароль
 
 
Рейтинг 4.68/50: Рейтинг темы: голосов - 50, средняя оценка - 4.68
АлександрКом
239 / 195 / 90
Регистрация: 21.10.2012
Сообщений: 904
1

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

29.03.2014, 23:39. Просмотров 9343. Ответов 20
Метки нет (Все метки)

Здравствуйте. Задача такая: есть массив данных (запись сигнала). В этом сигнале присутствуют 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, в данном случае есть ЦЕНТРАЛЬНЫЙ отсчет.

Скажите, будет ли работать данный алгоритм? Может нужно что исправить? Или есть ещё более простой способ?
0
Лучшие ответы (1)
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
29.03.2014, 23:39
Ответы с готовыми решениями:

Реализация ФНЧ Баттерворта
есть ли у кого уже готовый код расчета фильтра?

Расчет коэффициентов ФНЧ
как рассчитать входные коэффициенты для функции единичного отсчета и синусоида...

Результат моделирования ФНЧ
Ребята, имеется результат моделирования двух фильтров НЧ 3 и 5 порядков. Не...

ФНЧ при моделировании OFDM
Здравствуйте, Я моделирую OFDM-сигнал и при его приеме необходимо использовать...

Через что лучше ДПФ или ФНЧ?
короче есть задача сигнал пропускать через определенную частоту. То есть...

20
raxp
10188 / 6571 / 492
Регистрация: 28.12.2010
Сообщений: 21,166
Записей в блоге: 1
30.03.2014, 10:40 2
...через накопление и усреднение сигнала тоже можно.
0
vital792
2002 / 1274 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
30.03.2014, 10:49 3
Цитата Сообщение от АлександрКом Посмотреть сообщение
Преобразование Фурье входного сигнала -> преобразование Фурье коэффициентов КИХ-фильтра -> быстрая свертка КИХ-фильтра и сигнала -> обратное преобразование Фурье результатов свертки
Это все называется алгоритмом быстрой свертки. То что у вас в середине "быстрая свертка КИХ-фильтра и сигнала" заменяется просто умножением. Вот только при небольшом порядке фильтра, непосредственное выполнение свертки все таки выполняется быстрее этой цепочки. Легко даже посчитать, в каких случаях какой способ быстрее. По поводу вашего фильтра: преобразование фурье вашего фильтра - прямоугольник (если n бесконечно), так что его можно не выполнять, а просто обнулить нужный вам участок. Это простейший ких фильтр с прямоугольным окном, не самый лучший разумеется. Это не раз обсуждалось в данной ветке форума. Из за явления гиббса на границах полосы появляются пульсации. С ними борятся различными способами, самый простой из которых - домножение на весовое окно. Так что советую ничего не изобретать, а просто использовать обычный фильтр и обычную свертку. пример
0
eSergey76
0 / 0 / 0
Регистрация: 26.12.2012
Сообщений: 5
30.03.2014, 15:06 4
Более простой способ - WinFilter08.zip. Выбираете тип фильтра и генерите С или VHDL код. Все.
0
АлександрКом
239 / 195 / 90
Регистрация: 21.10.2012
Сообщений: 904
30.03.2014, 16:34  [ТС] 5
Цитата Сообщение от vital792 Посмотреть сообщение
Так что советую ничего не изобретать, а просто использовать обычный фильтр и обычную свертку. пример
Посмотрел пример по ссылке, там вы приводите свой код. Хотел спросить: в коде, я так понял, создается ИХ фильтра и происходит свертка её с сигналом. На выходе получаем отфильтрованный сигнал,да? Хотел узнать по поводу некоторых переменных, правильно ли я понял, что Wn - это частота дискретизации, order - порядок фильтра? Только не понял, наверное, самый главный момент - т.к. это ФНЧ, то должна задаваться частота среза, а в коде как-то не "осознал" этого момента.. Не могли бы пояснить?
0
vital792
2002 / 1274 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
30.03.2014, 17:37 6
Wn - как раз частота среза. Частота дискретизации тут не важна. Точнее она не используется, т.к. тут нормированная частота, то есть частоту среза задаете от 0 до 1, считая за 1 половину частоты дискретизации
0
АлександрКом
239 / 195 / 90
Регистрация: 21.10.2012
Сообщений: 904
30.03.2014, 20:00  [ТС] 7
Попробовал ваш фильтр при демодуляции ОБП сигнала, что- то не так видно задаю, т.к. на выходе фильтра одни нули. (с 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;
}
0
vital792
2002 / 1274 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
30.03.2014, 20:41 8
Лучший ответ Сообщение было отмечено raxp как решение

Решение

Цитата Сообщение от АлександрКом Посмотреть сообщение
(с 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;
}
2
АлександрКом
239 / 195 / 90
Регистрация: 21.10.2012
Сообщений: 904
30.03.2014, 22:44  [ТС] 9
Спасибо, теперь стало понятнее Остался последний вопрос: если судить по книжке, по которой я делал свой сигнал, там после демодуляции и отфильтровки сигнал получается такой же, как и тот, который я использовал для модуляции (s_M), только по амплитуде в 2 раза меньше. У меня же в коде, при сравнении значений s_M и result, такой зависимости нет.. В чем может быть проблема?
0
vital792
2002 / 1274 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
30.03.2014, 23:18 10
Цитата Сообщение от АлександрКом Посмотреть сообщение
В чем может быть проблема?
в крутизне спада характеристики фильтра. Для такого маленького порядка она довольно плавная и вч колебания фильтруются не полностью.
Вот ваш исходный сигнал:
Реализация ФНЧ

А вот результат после фильтрации:
Реализация ФНЧ

Если порядок фильтра увеличить до 40, ачх будет круче и результат будет такой:
Реализация ФНЧ
2
АлександрКом
239 / 195 / 90
Регистрация: 21.10.2012
Сообщений: 904
01.04.2014, 15:23  [ТС] 11
Спасибо, за объяснение, теперь, действительно стало понятно

Добавлено через 2 часа 7 минут
И снова неприятность я думал, что порядок фильтра задается #define order 20, но при замене 20 на 40 получил практически тот же результат, что и у вас, но до замены порядка фильтра (когда вч колебания фильтруются не полностью)
0
vital792
2002 / 1274 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
01.04.2014, 17:37 12
да верно, у меня вышла ошибочка. Вместо
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 умножается все выражение.
2
АлександрКом
239 / 195 / 90
Регистрация: 21.10.2012
Сообщений: 904
02.04.2014, 20:04  [ТС] 13
Спасибо большое
0
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
02.05.2015, 15:08 14
Реализация полосового фильтра Баттерворта 3 порядка:
ФНЧ-прототип:
http://www.cyberforum.ru/cgi-bin/latex.cgi?H(S)=\frac{\beta ^{3}}{S^{3}+2\beta S^2+2\beta ^{2}S+\beta^3}
где
http://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-преобразования:
http://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}

где
http://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 - частота Найквиста.

http://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 минут
Если нужен обычный ФНЧ, то расчёты существенно упрощаются:
ФНЧ-прототип:
http://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 порядка:
http://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}
http://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
raxp
10188 / 6571 / 492
Регистрация: 28.12.2010
Сообщений: 21,166
Записей в блоге: 1
02.05.2015, 20:46 15
...Матлабом не пользуюсь.

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

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

Я конечно, мог бы формулы окончательные для коэффициентов фильтра
уже есть онлайн-генератор готового Си-кода http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html
0
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
04.05.2015, 11:22 16
ПФ 4-го порядка.
Матрица z-преобразования:
http://www.cyberforum.ru/cgi-bin/latex.cgi?<br />
\begin{bmatrix}<br />
1 & 1 &1  & 1 & 1 & 0 & 0 &0  & 0\\ <br />
 4\Psi & 3\Psi  &  2\Psi &\Psi   &0  & 0 & 0 & 0 &0 \\ <br />
 4&  2& 0 & -2 & -4 & 0 & \Psi^2 &3\Psi^2   &6\Psi^2 \\ <br />
 12\Psi & 3 & -2\Psi  &-3\Psi   &0  & 0 &0  &\Psi^3   & 4\Psi^3 \\ <br />
 12\Psi^2 &0  &-2  & 0 &6     &0  &-2\Psi^2  &0 &\Psi^4+6 \\ <br />
12\Psi & -3 & -2\Psi  &3\Psi   &0  & 0 &0  &-\Psi^3   & 4\Psi^3 \\<br />
 4&  -2& 0 & 2 & -4 & 0 & \Psi^2 &-3\Psi^2   &6\Psi^2 \\ <br />
  4\Psi & -3\Psi  &2\Psi   & - \Psi &  0&  0& 0 &  0&0 \\ <br />
 1&  -1& 1 &  -1& 1 & 0 & 0 & 0 &0 <br />
\end{bmatrix}\begin{bmatrix}<br />
1\\ <br />
v\beta \\ <br />
w\beta^2 \\ <br />
v\beta^3 \\ <br />
\beta^4 \\ <br />
v\beta^3 \\ <br />
w\beta^2 \\ <br />
v\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 />
A_7\\ <br />
A_8<br />
\end{bmatrix}
Обозначения, как в предыдущем посте, недостающие коэффициенты в статье http://www.dsplib.ru/content/filters/butterex/butterex.html

Добавлено через 17 часов 5 минут
Матрица z-преобразования для фильтра 15 порядка вот так выглядит:
1111111111111111
-15-13-11-9-7-5-3-113579111315
105775333175-3-7-7-3517335377105
-455-273-143-57-715177-7-17-15757143273455
136563722121-43-35-32121-3-35-43212216371365
-3003-1001-1439977-1-39-2121391-77-9914310013003
50051001-143-187-116525-35-352565-11-187-14310015005
-6435-42942999-99-454535-35-454599-99-4294296435
6435-429-4299999-45-453535-45-459999-429-4296435
-50051001143-1871165-25-353525-65-11187-143-10015005
3003-100114399-77-139-21-2139-1-7799143-10013003
-1365637-2212143-35321-21-335-43-21221-6371365
455-273143-57715-1777-17157-57143-273455
-10577-5333-1753-77-3-517-3353-77105
15-1311-97-53-1-13-57-911-1315
-11-11-11-11-11-11-11-11
               
Столбцы матрицы - это коэффициенты произведения полиномов вида (z-1)^n(z+1)^m , где n,m=0...N, n+m=N (порядок фильтра)

Добавлено через 4 часа 46 минут
Внимание!
В матрице z-преобразования ПФ 3 порядка опечатка - изменить знак каждого элемета последней строки.
В матрице z-преобразования ПФ 4 порядка опечатка - 2 столбец, строки 4 и 6: Значения http://www.cyberforum.ru/cgi-bin/latex.cgi?3\Psi ;  -3\Psi.
1
Red_meat
0 / 0 / 0
Регистрация: 14.05.2015
Сообщений: 9
14.05.2015, 20:10 17
Добрый человек, а можно к Вам с еще одним вопросом?
Почему когда Wn = 0.05 это ФНЧ, Wn = 0.5 ФВЧ, а при большинстве прочих выражений это полосной фильтр? И почему это работает, а то что написано у "Айфичер Э. - Цифровая обработка сигналов. Практический подход. 2-е издание. 2004" в методе взвешивания нет?

Добавлено через 21 минуту
Если быть точнее, мне нужно регулировать частоту среза и ширину полосы перехода. Может подскажите как это сделать или что я делаю не так, когда сворачиваю сигнал с фильтром?
0
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
14.05.2015, 22:25 18
Я если честно - вопроса не понял совсем
В коэффициентах (для фильтра Баттерворта) очень легко переход ФНЧ-ФВЧ получить. Надо создать ФНЧ с частотой (1-Fs) ( в относительных частотах) и потом поменять знаки у нечётных коэффициентов. С полосовыми фильтрами - сложнее. Могу посоветовать ссылку на работы Костантинидеса.
Приведите конкретный пример фильтра, чтобы я понял суть вопроса
1
Red_meat
0 / 0 / 0
Регистрация: 14.05.2015
Сообщений: 9
14.05.2015, 23:22 19
Я не про Баттерворта. Я про начала темы и код в ней. Там написан ФНЧ с заранее указанной частотой среза, порядком и прочим. Если хотите - посмотрите.
Что пытаюсь сделать я, генерировать частотные фильтры по типу(ФНЧ, ФВЧ, полосной, заграждающий) частоте среза и затуханию дБ в полосе. Книжка "Айфичер Э. - Цифровая обработка сигналов." говорит о методе взвешивания.
Гласит он, смешать две таблицы(прилагаются) линейной сверткой, и будет вам счастье. Порядок фильтра предлагают брать нечетный, а коэффициенты симметричные. Вроде сказка, а не метод, а получается фигня.

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с, и по частотам будто узким окном проходят. При одной это ФНЧ, при другой ФВЧ, а при всех остальных полосной. Но где он окажется, этот не задавленный фильтром участок частот, я не знаю.

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

Ну и если уж пошла такая пьянка, может кто подскажет, как АЧХ такого фильтра посмотреть в матлабе можно?
0
Миниатюры
Реализация ФНЧ   Реализация ФНЧ  
Red_meat
0 / 0 / 0
Регистрация: 14.05.2015
Сообщений: 9
16.05.2015, 02:34 20
Задача решилась! Оставлю код, может кому поможет.
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;
}
0
16.05.2015, 02:34
MoreAnswers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
16.05.2015, 02:34

Расчет порядка нормированного прототипа ФНЧ при расчете полосового фильтра
Добрый день! После длительного перерыва опять занялся разработкой приложения...

ФНЧ третьего порядка, коэффициент передачи в полосе пропускания 6 дБ.
Доброго времени суток. Требуется разработать схему ФНЧ третьего порядка, одним...

? реализация сервопривода
Здравствуйте, так и не нашел как связаться с di-halt, напишу тут. Как воздух...


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

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

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Рейтинг@Mail.ru