1차 lpf 를 설계하고 구현해 보겠다.
1. 개념
신호의 낮은 주파수 부분을 통과시키는 것이다.
먼저 1차 lpf 의 전달 함수는 다음과 같이 정의 할 수 있다.
$$ H(s)=\frac{1}{1+RCs} $$
이 수식이 정말로 낮은 주파수 부분은 통과시키고 높은 주파수 부분은 통과시키지 않는 것은 다음과 같은 matlab 코드로 확인할 수 있다.
% 시간 상수 tau 정의
tau = 1; % 여기서는 tau = 1초로 설정
% 전달함수 정의
s = tf('s');
H = 1 / (tau * s + 1);
% Bode plot 그리기
figure;
bode(H);
grid on;
title('1차 저주파 필터의 Bode Plot');
결과에서 볼 수 있듯이 높은 주파수로 갈 수록 이득이 떨어지는 것을 볼 수 있다.
2. 수식 유도
전달 함수
$$ H(s)=\frac{1}{1+RCs} $$
시정수 $\tau$ 는 RC 로 정의되므로 다시 정리하면
$$ H(s) = \frac{1}{\tau s+1} $$
또 다른 표현으로는
$$ H(s) = \frac{\omega_c}{s + \omega_c} $$
$$ \omega_c = \frac{1}{\tau} $$
코드로 변환하기 위해 H(s) 에 Bilinear Transform 적용
$$
s = \frac{2}{T} \cdot \frac{1 – z^{-1}}{1 + z^{-1}}
$$
Bilinear Transform을 적용하면 s-도메인 전달 함수를 z-도메인으로 변환할 수 있습니다. 이를 위해 s를 다음과 같이 대체합니다:
여기서 T는 샘플링 주기입니다. 전달 함수 H(s) 를 변환하면 다음과 같이 됩니다: $$ H(z) = H\left( \frac{2}{T} \cdot \frac{1 – z^{-1}}{1 + z^{-1}} \right) $$ 구체적인 변환 과정: $$ H(z) = \frac{\omega_c}{\frac{2}{T} \cdot \frac{1 – z^{-1}}{1 + z^{-1}} + \omega_c} $$ $$ H(z) = \frac{\omega_c (1 + z^{-1})}{\frac{2}{T} (1 – z^{-1}) + \omega_c (1 + z^{-1})} $$ $$ H(z) = \frac{\omega_c (1 + z^{-1})}{\left(\frac{2}{T}\right) (1 – z^{-1}) + \omega_c (1 + z^{-1})} $$ $$ H(z) = \frac{\omega_c (1 + z^{-1})}{\left(\frac{2}{T}\right) – \left(\frac{2}{T}\right) z^{-1} + \omega_c + \omega_c z^{-1}} $$ $$ H(z) = \frac{\omega_c (1 + z^{-1})}{\left(\frac{2}{T} + \omega_c\right) + \left(\omega_c – \frac{2}{T}\right) z^{-1}} $$ 여기서, $a = \frac{2}{T} + \omega_c$ 와 $b = \omega_c – \frac{2}{T}$ 로 정의하면:
$$ H(z) = \frac{\omega_c (1 + z^{-1})}{a + b z^{-1}} $$
or
$$
H(s) = \frac{1}{\tau s + 1}
$$
$$
H(z) = H\left( \frac{2}{T} \cdot \frac{1 – z^{-1}}{1 + z^{-1}} \right)
$$
변환 과정을 자세히 살펴보면:
$$
H(z) = \frac{1}{\tau \left( \frac{2}{T} \cdot \frac{1 – z^{-1}}{1 + z^{-1}} \right) + 1}
$$
$$
H(z) = \frac{1}{\frac{2\tau}{T} \cdot \frac{1 – z^{-1}}{1 + z^{-1}} + 1}
$$
분모를 정리하면:
$$
H(z) = \frac{1}{\frac{2\tau (1 – z^{-1})}{T (1 + z^{-1})} + 1}
$$
$$
H(z) = \frac{T (1 + z^{-1})}{2\tau (1 – z^{-1}) + T (1 + z^{-1})}
$$
$$
H(z) = \frac{T (1 + z^{-1})}{T + 2\tau – 2\tau z^{-1} + T z^{-1}}
$$
여기서 $a = T + 2\tau$이고 $b = T – 2\tau$라고 정의하면:
$$
H(z) = \frac{T (1 + z^{-1})}{a + b z^{-1}}
$$
위의 전달 함수를 시간 영역으로 변환하면 다음과 같은 차분 방정식을 얻을 수 있습니다:
$$
H(z) = \frac{Y(z)}{X(z)} = \frac{T (1 + z^{-1})}{a + b z^{-1}}
$$
이를 시간 영역으로 풀면:
$$
Y(z) \left(a + b z^{-1}\right) = X(z) \left(T (1 + z^{-1})\right)
$$
여기서 $Y(z)$와 $X(z)$는 각각 출력과 입력의 z-변환입니다.
시간 영역으로 변환하면:
$$
y[n] a + y[n-1] b = x[n] T + x[n-1] T
$$
이를 $y[n]$에 대해 풀면:
$$
y[n] = \frac{T x[n] + T x[n-1] – b y[n-1]}{a}
$$
3. 코드 구현
위에 식을 c++ 코드로 변환하면
만약 아래 식으로 시작했다면 $$ H(s) = \frac{1}{\tau s+1} $$
class LPF {
public:
LPF(float sampleRate, float cutoffFrequency) {
T = 1.0 / sampleRate;
omega_c = 2 * M_PI * cutoffFrequency;
tau = 1.0 / omega_c;
a = T + 2 * tau;
b = T - 2 * tau;
prevInput = 0.0;
prevOutput = 0.0;
}
float process(float input) {
float output = (T * (input + prevInput) - b * prevOutput) / a;
prevInput = input;
prevOutput = output;
return output;
}
private:
float T;
float omega_c;
float tau;
float a, b;
float prevInput;
float prevOutput;
};or
만약 아래 식으로 시작했다면
$$ H(s) = \frac{\omega_c}{s + \omega_c} $$
class LPF {
public:
LPF(float sampleRate, float cutoffFrequency) {
T = 1.0 / sampleRate;
omega_c = 2 * M_PI * cutoffFrequency;
a = (2.0 / T) + omega_c;
b = omega_c - (2.0 / T);
prevInput = 0.0;
prevOutput = 0.0;
}
float process(float input) {
float output = (omega_c * (input + prevInput) - b * prevOutput) / a;
prevInput = input;
prevOutput = output;
return output;
}
private:
float T;
float omega_c;
float a, b;
float prevInput;
float prevOutput;
};4. 최종 구현 결과 확인
sample data 을 생성해서 이를 filter 적용한 후 sample data와 filter 적용된 data 각각 txt로 저장하는 코드이다.
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <string>
class LPF {
public:
LPF(float sampleRate, float cutoffFrequency) {
T = 1.0 / sampleRate;
omega_c = 2 * M_PI * cutoffFrequency;
tau = 1.0 / omega_c;
a = T + 2 * tau;
b = T - 2 * tau;
prevInput = 0.0;
prevOutput = 0.0;
}
float process(float input) {
float output = (T * (input + prevInput) - b * prevOutput) / a;
prevInput = input;
prevOutput = output;
return output;
}
private:
float T;
float omega_c;
float tau;
float a, b;
float prevInput;
float prevOutput;
};
std::vector<double> filterData(const std::vector<double>& data, float cutoffFrequency, float sampleRate, size_t startIdx, size_t endIdx) {
LPF lpf(sampleRate, cutoffFrequency);
std::vector<double> filteredData(data.size(), 0.0);
for (size_t step = startIdx; step <= endIdx; ++step) {
filteredData[step] = lpf.process(data[step]);
}
return filteredData;
}
void saveData(const std::vector<double>& data, const std::string& filePath) {
std::ofstream outputFile(filePath);
if (!outputFile.is_open()) {
std::cerr << "파일을 열 수 없습니다!" << std::endl;
return;
}
for (const auto& value : data) {
outputFile << value << "\n";
}
outputFile.close();
}
int main() {
float sampleRate = 250.0; // Sample rate in Hz
float cutoffFrequency = 3.0; // Cutoff frequency in Hz
// 큰 주파수와 작은 주파수를 더한 샘플 데이터 생성
std::vector<double> data(6000);
for (size_t i = 0; i < data.size(); ++i) {
double largeSine = sin(2 * M_PI * 0.5 * i / sampleRate); // 0.5 Hz의 사인파
double smallSine = 0.5 * sin(2 * M_PI * 5 * i / sampleRate); // 5 Hz의 사인파 (작은 진폭)
data[i] = largeSine + smallSine;
}
// 각 열에 대해 3000번부터 5000번까지의 데이터에 1차 저역 통과 필터 적용
size_t startIdx = 3000;
size_t endIdx = 5000;
size_t dataSize = data.size();
if (dataSize < endIdx) {
std::cerr << "데이터의 길이가 충분하지 않습니다!" << std::endl;
return 1;
}
std::vector<double> filteredData = filterData(data, cutoffFrequency, sampleRate, startIdx, endIdx);
// 샘플 데이터와 필터링된 데이터를 파일에 저장
saveData(data, "./sample_data.txt");
saveData(filteredData, "./filtered_data.txt");
return 0;
}
위에가 cpp 코드
% sample_data.txt와 filtered_data.txt 파일 경로
sampleDataFilePath = './sample_data.txt';
filteredDataFilePath = './filtered_data.txt';
% 데이터 로드
sampleData = load(sampleDataFilePath);
filteredData = load(filteredDataFilePath);
% 데이터 범위 설정 (3000번부터 5000번까지)
startIdx = 3000;
endIdx = 5000;
% 그래프 그리기
figure;
plot(sampleData(startIdx:endIdx), 'b-');
hold on;
plot(filteredData(startIdx:endIdx), 'r', 'LineWidth',5);
hold off;
% 그래프 제목 및 라벨 설정
title('Original and Filtered Data');
xlabel('Sample Index');
ylabel('Amplitude');
legend('Original Data', 'Filtered Data');
grid on;위에가 결과 확인을 위하 matlab 코드
위에 가 최종 결과이다.
추후 다른 필터도 구현해 보겠다
0 Comments