C++을 이용하여 LPF를 설계하는 과정을 했던 프로젝트이다. 예전에 했던프로젝트라 부족한면도 많긴 하지만, 코딩하면서 꽤 재미를 느끼게 된 프로젝트이고, 코딩과 신호처리를 한번에 아우를 수 있어서 뜻깊었던 프로젝트였던 것 같다. 살면서 처음해본 프로젝트이기도 하다.
언어는 C++를 사용했고, mfc나 기타 헤더파일없이 직접 한땀한땀 쳐가면서 해본 프로젝트라서 나름 감명깊은 프로젝트였다. 적분기호들로 나타내져있는 수식들을 디지털로 나열하고, 배열의 공간 개념등이 어느정도 적용되어 있어야 하며, Foruier Transform 및 기타 신호처리 지식이 어느정도 있으면 디지털에서의 신호필터설계가 어떻게 되는지 이해하기 굉장히 편할 것이라고 생각한다.
우선 LPF란 무엇인지 아는것이 필요하다.
LPF란 Low Pass Filter 즉, 낮은 주파수 대역을 모두 날려버리는 필터이다. 이를 어떤식으로 처리하는지 간단히 알아보자.
#define SizeofPulse 512 // 총 관찰하는 구간
#define N_pulse 512 // N으로 잡은 값
#define SizeofSINC 512 //생성할 SINC 범위
#define pi 3.1415
void clear_filter(double pulse[]); // pulse 소규모 작은값을 정리해주는 필터
void Fourier_Transform(double Pulse[],double Yn[]); // Pulse를 F.T하여 Yn에 저장
void convolution(double Xn[], double Hn[], int length,double Yn[]); // Xn, Yn을 convolution하여 Yn에 저장
void S_pulse(double SIN[], int Hz); // Hz를 주기로 가지는 sin파형 생성
void Sum_pulse(double pulse1[], double pulse2[], double pulse3[]); // pulse1,pulse2를 더한값을 pulse3에 저장
void R_pulse(int half_width, int period, double Pulse[]); // 반복 구형파 주기함수 생성
void Out_pulse(double pulse[]); //pulse 출력함수
void SINC_pulse(double pulse[],double cut_fre); // SINC파형 생성 (cut_fre : 차단주파수)
void Butterworth_pulse(double pulse[], double cut_fre, int n); // cut_fre : 차단주파수, n : 곱해주는 차수
void Chebyshev_pulse(double pulse[], int N, double delta,double cut_fre); // N : 스윙횟수 , delta : Ripple 폭, cut_fre : 차단주파수
void Whitenoise_pulse(double pulse[], int range);생성한 함수들과 정의들은 다음과 같다. 총 512개의 샘플을 관찰하였고, 각 펄스를 생성하는 함수들을 생성하였다.
각 함수들은 아래와 같이 만들어진다. 크게 보면, for문을 사용하여 512개의 배열을 1초라고 본 후에, 그안에 진폭, 주파수, 그리고 sin인지 cos인지, sinc인지 등등에 따라 수식을 넣어주었다고 생각하면 된다.
보통 신호에서 x(t) = sin ( 2 pi f t) 의 형식으로 많이 적기에 그대로 반영해보도록 노력했다.
void Whitenoise_pulse(double pulse[], int range)
{
for (int time = 0; time < SizeofPulse; time++)
{
pulse[time] = ((rand() % range) - range / 2);
}
}
void Chebyshev_pulse(double pulse[], int N, double delta, double cut_fre)
{ // N : 스윙횟수 , delta : Ripple 폭, cut_fre : 차단주파수
double square_epsilon = 1 / delta - 1;
double Cn;
for (int k = 0; k < SizeofSINC; k++) {
if (k > cut_fre)
Cn = cosh(N * acosh(k / cut_fre));
else
Cn = cos(N * acos(k / cut_fre));
pulse[k] = 1 / (sqrt(1 + square_epsilon * Cn * Cn));
}
}
void Butterworth_pulse(double pulse[], double cut_fre, int n)
{ // cut_fre : 차단주파수, n : 곱해주는 차수
for (int k = 0; k < SizeofSINC; k++) {
pulse[k] = 1 / (sqrt(1 + pow((k / cut_fre), n)));
}
}
void SINC_pulse(double pulse[],double cut_fre) //cut_fre를 차단주파수로 하는 sinc함수 생성
{
pulse[0] = 1;
for (int k = 1; k < SizeofSINC; k++) {
pulse[k] = (sin(( k *2*cut_fre) / SizeofSINC) )/ ( k * 2 * cut_fre / SizeofSINC);
}
for (int k = SizeofSINC; k < SizeofPulse; k++) {
pulse[k] = 0;
}
}
void clear_filter(double pulse[])
{
for (int i = 0; i < SizeofPulse; i++) {
if ((pulse[i] < 5) && (pulse[i] > -5)) {
pulse[i] = 0;
}
}
}
void Fourier_Transform(double Pulse[],double Yn[])
{
int time = 0;
double Re_Xf[SizeofPulse]; //실수부분
double Im_Xf[SizeofPulse]; //허수부분
for (int k = 0; k < SizeofPulse; k++)
{
double sum_Re = 0, sum_Im = 0;
for (time = 0; time < SizeofPulse; time++)
{
sum_Re += Pulse[time] * cos((pi * 2 * k * time) / N_pulse);
sum_Im += Pulse[time] * sin((pi * 2 * k * time) / N_pulse);
}
Re_Xf[k] = sum_Re;
Im_Xf[k] = -sum_Im;
}
for (int k = 0; k < SizeofPulse; k++)
{
clear_filter(Im_Xf);
Yn[k]= sqrt( (Re_Xf[k])* (Re_Xf[k]) + (Im_Xf[k]) * (Im_Xf[k])) ; //크기 출력
// cout << -atan((Im_Xf[k] / Re_Xf[k])) * (180 / pi) << "\n"; //phase 출력 (각도)
}
}
void convolution(double Xn[], double Hn[],int length, double Yn[])
{ // Xn, Yn을 convolution하여 Yn에 저장
int time = 0;
double yn[2* SizeofPulse] = { };
for (int k = 0; k < length+length; k++)
{
yn[k] = 0;
for (time = 0; time < length+length; time++)
{
if ((k - time) >= 0)
{
if( (k - time) < length && time < length)
{
yn[k] += Xn[time] * Hn[k - time];
}
}
}
}
for (int k = 0; k < length; k++) {
Yn[k] = yn[k];
}
}
void S_pulse(double SIN[], int Hz)
{
for (int i = 0; i < SizeofPulse; i++) {
SIN[i] = sin(((2 * pi * (i * Hz + 0))) / N_pulse);
} //0에 숫자를 넣으면 그만큼 위상변화
}
void Sum_pulse(double pulse1[], double pulse2[], double pulse3[]) //두개의 펄스를 더해주는 함수
{
for (int i = 0; i < SizeofPulse; i++) {
pulse3[i] = pulse1[i] + pulse2[i];
}
}
void R_pulse(int half_width, int period, double Pulse[]) // 주기 사각펄스 생성
{
for (int time = 0; time < SizeofPulse; time++)
{
if (time <= half_width)
Pulse[time] = 1;
else if (time % period <= half_width)
Pulse[time] = 1;
else if (period - (time % period) <= half_width)
Pulse[time] = 1;
else
Pulse[time] = 0;
}
}
void Out_pulse(double pulse[]) //출력함수
{
for (int i = 0; i < SizeofPulse; i++) {
cout << pulse[i] << "\n";
}
}
chebyshev 필터와 butterworth필터가 나오는데, 이는 조금 후에 설명하도록 하겠다.
푸리에에 대해선 많이 다루지 않겠지만, 한 시스템을 지나간 결과가 y(n) = h(n) * x(n) 이라는 사실은 알고 있을것이다. 여기서 *은 곱이 아닌, 컨볼루션이다. 이를 각자 주파수축에서 바꾸면, Y(f) = H(f) * X(f) 로 바뀌게 된다. 여기서의 *은 단순 곱이다. 시스템분석을위해서 컨볼루션하는것도 해보고, 주파수측으로 바꾼후 단순곱도 해보면서 여러가지 실험을 해볼 수있을 것이다.
void main()
{
double SIN_30[SizeofPulse]; // sin(1Hz)
double SIN_40[SizeofPulse]; // sin(50Hz)
double SIN_50[SizeofPulse];
double SIN_60[SizeofPulse];// sin(100Hz)
double Sum_SIN[SizeofPulse]; // sin (1Hz) + sin(50Hz)
S_pulse(SIN_30, 30);
S_pulse(SIN_40, 40);
S_pulse(SIN_50, 50);
S_pulse(SIN_60, 60);
double ex1[SizeofPulse] = { 1,1,1,1,1, };
double SIN[SizeofPulse];
for (int i = 0; i < SizeofPulse; i++)
{
SIN[i] = SIN_30[i] + SIN_40[i] + SIN_50[i] + SIN_60[i];
}
Sum_pulse(SIN_30, SIN_40, Sum_SIN);
double R_Pulse[SizeofPulse]; // 사각주기파형 생성
R_pulse(30, 200, R_Pulse);
double SINC[SizeofPulse]; //SINC 생성
double Butterworth[SizeofPulse]; // Butterworth 필터
double Chebyshev[SizeofPulse]; //Chebyshev 필터
SINC_pulse(SINC, 50 * (double(SizeofSINC) / double(N_pulse))); //SINC를 다 안가져가도 차단주파수가 바뀌지 않게함
Butterworth_pulse(Butterworth, 50, 10);
Chebyshev_pulse(Chebyshev, 100, 0.9, 50);
double Whitenoise[SizeofPulse];
double Sum_SIN_Noise[SizeofPulse];
Whitenoise_pulse(Whitenoise, 5);
Sum_pulse(Sum_SIN, Whitenoise, Sum_SIN_Noise);
double Yn[SizeofPulse];
double Hn[SizeofPulse];
Fourier_Transform(Sum_SIN_Noise,Hn); //Hn = F{ H(f) }
convolution(SINC, SIN,SizeofPulse,Yn); //Yn = Xn convolution Hn
// Fourier_Transform(Yn,Yn); //Y(f) = F{Y(n)}
Out_pulse(SINC);
}다음과 같은 main문을 작성할 수 있었다. 쉽게생각하면, 여러 함수들을 다 만들어두고, Yn,Hn공간만 놔둔 ㅎ후에, 푸리에를 진행하고 컨볼루션도 하고, 주파수에서 단순곱한다음에 지지고 복고 하면서 최종적으로 나오게 되는 Out_pulse를 사용하여 그 값을 자유자제로 찍을 수 있는 프로그램이다. mfc를 적용한 프로그램이 아니라 보는데 불편하고, 수치를 그래프로 그리려면 직접 다른 프로그램을 이용해야하는 불편함이 있긴 하다.
다시 신호처리적 내용으로 돌아와서, 먼저 butterworth필터를 보겠다.
위 신호가 butterworth를 Fourier Transfom하여 나타낸 형태이다.
위가 시간축 측면이며, 디지털 필터특성상 음수값이 절대값으로 치환되어 나오게 설정하였다. sinc와 굉장히 유사하다고도 생각할 수있다.
위 신호를 sin1Hz에 sin 50Hz 두 신호가 더해진, 값을 넣고, 이 butterworth 필터를 H(n)인 시스템에 집어넣으면,
위와 같은 결과를 얻을 수 있다. 보면 알겠지만, sin 1Hz성분은 거의 그대로 남아있고, 50Hz성분이 완전히 제거되진 않았지만, 많이 제거되어 있음을 알 수 있다.
위에서 보면 더 알기 편할 것이다. 50Hz부분의 주파수의 진폭값이 깎여있는것을 알 수있다. 원래신호는 둘의 크기가 거의 비슷하다고 생각하면 된다. 사실 주파수측면에서는 단순 곱이기 때문에 머릿속으로 조금만 생각해도 바로 답이 나올 것이다.
위와 비슷한 필터로 chbyshev필터가 있는데 이 필터는 아래와 같이 생겼다.
굉장히 rectangle신호와 비슷하나, 위에서 스윙되는 성질이 있음을 알 수있다.
위 필터 역시 sinc나 butterworth필터와 굉장히 유사하게 생겼다.
이역시 결과는 50Hz성분이 날라갔음을 확인할 수있다.
그렇다면 이런 질문이 하나 들 수있을것이다 그냥 주파수에서의 Rectangle신호, 즉 sinc파형을 LPF필터로 사용하면 안되는가?
당연히 들 수있는 의문이고, 프로그램상으로도 간단하기에 구현해보았다.
역시 잘됨을 알 수있다. 하지만, 위에서 설명하였듯이 이는 그저 rectangle신호를 구현하여 단순곱을 한 후, 다시 I.Fourier를 통하여 복구한것이기에 잘된거고, sinc파형을 그대로 convolution한다면,아래와 같이 될 것이다.
우선, sinc를 800개와 5개 샘플링을 해가면서 생각을 해보았다.(참고로 위의 필터들은 모두 512개를 기준으로 하였다.)
결과를 보면, sinc파형이 오히려 필터효과가 더 덜 된다는 것을 알 수있다.
'프로젝트 > 소규모프로젝트들' 카테고리의 다른 글
| JS_리액트_TodoList (0) | 2021.12.27 |
|---|---|
| MFC를 이용한 디지털신호처리 프로그램 (0) | 2021.12.15 |
| JS_독서관리 시스템 (2) | 2021.12.10 |
| GitHub 클론코딩-2 (0) | 2021.11.27 |
| GitHub 클론코딩-1 (0) | 2021.11.26 |