1. EWT의 이론적 기초와 핵심 개념

1.1 기존 신호 분해 기법과의 비교

1.2 EWT의 주요 혁신 요소

• 주파수 도메인 자동 분할: 푸리에 스펙트럼의 국소 최대점을 기준으로 주파수 대역을 동적으로 나눔
• 직교 웨이블릿 필터 설계: 각 대역에 대해 Meyer 웨이블릿 기반 필터 생성
• 무손실 재구성 보장: 원본 신호를 완벽하게 복원 가능
• 수학적 해석 가능성: EMD와 달리 수학적으로 엄밀히 정의됨

2. MATLAB에서의 구현

2.1 EWT 메인 함수 (empirical_wavelet_transform.m)

function [components, boundaries] = empirical_wavelet_transform(input_signal, varargin)
% 경험적 웨이블릿 변환
% 입력:
%   input_signal - 1차원 신호 벡터
%   varargin     - 'numPeaks', 'logScale', 'boundaryFactor', 'waveletWidth'
% 출력:
%   components   - 추출된 구성 성분들 (행 단위 저장)
%   boundaries   - 주파수 분할 경계

% 기본 파라미터 설정
numPeaks = 4;
useLogScale = true;
boundaryFactor = 0.4;
waveletWidth = 0.1;

% 가변 인자 처리
for i = 1:2:length(varargin)
    switch lower(varargin{i})
        case 'numpeaks', numPeaks = varargin{i+1};
        case 'logscale', useLogScale = varargin{i+1};
        case 'boundaryfactor', boundaryFactor = varargin{i+1};
        case 'waveletwidth', waveletWidth = varargin{i+1};
    end
end

signalLength = length(input_signal);
halfLength = ceil(signalLength / 2);

% 푸리에 변환 수행
spectrum = fft(input_signal);
magnitude = abs(spectrum(1:halfLength));

% 로그 스케일 적용 시 해상도 향상
if useLogScale
    processedMagnitude = log(magnitude + 1e-10);
else
    processedMagnitude = magnitude;
end

% 주파수 축 계산
frequencyAxis = (0:halfLength-1) * (1/signalLength);

% 주파수 대역 자동 분할
boundaries = detect_frequency_boundaries(processedMagnitude, frequencyAxis, numPeaks, boundaryFactor);

% 경험적 웨이블릿 필터 생성
filterBank = construct_filter_bank(boundaries, signalLength, waveletWidth);

% 필터링을 통한 성분 추출
numComponents = length(boundaries) - 1;
components = zeros(numComponents, signalLength);

for idx = 1:numComponents
    filteredSpectrum = spectrum .* filterBank(idx, :)';
    components(idx, :) = real(ifft(filteredSpectrum));
end
end

%% 주파수 경계 탐지 함수
function bounds = detect_frequency_boundaries(mag, freq, peakCount, factor)
L = length(mag);
[peaks, locations] = findpeaks(mag, 'MinPeakDistance', round(L/20));

% 주요 피크만 선택
if length(peaks) > peakCount
    [~, sortedIdx] = sort(peaks, 'descend');
    selectedLocs = locations(sortedIdx(1:peakCount));
else
    selectedLocs = locations;
end

[sortedLocs, ~] = sort(selectedLocs);
bounds = zeros(length(sortedLocs)+2, 1);
bounds(1) = 1;

for k = 1:length(sortedLocs)-1
    midPoint = (sortedLocs(k) + sortedLocs(k+1)) / 2;
    delta = factor * (sortedLocs(k+1) - sortedLocs(k)) / 2;
    candidate = sortedLocs(k) + delta;
    bounds(k+1) = max(candidate, bounds(k) + 1);
end
bounds(end) = L;
bounds = unique(round(bounds));
end

%% 필터 뱅크 구성 함수
function bank = construct_filter_bank(edges, N, width)
numBands = length(edges) - 1;
bank = zeros(numBands, N);
freqAxis = (0:N-1) / N;

for b = 1:numBands
    w_low = (edges(b) - 1) / N;
    w_high = (edges(b+1) - 1) / N;
    
    for n = 1:N
        w = freqAxis(n);
        if w <= w_low
            bank(b, n) = 0;
        elseif w >= w_high
            bank(b, n) = 0;
        else
            normalized = (w - w_low) / (w_high - w_low);
            bank(b, n) = cos(pi/2 * transition_function(width * normalized));
        end
    end
    
    % 복소수 대칭성 유지
    bank(b, N/2+2:end) = fliplr(bank(b, 2:N/2));
end

% 에너지 정규화
for b = 1:numBands
    maxVal = max(abs(bank(b, :)));
    if maxVal > 0
        bank(b, :) = bank(b, :) / maxVal;
    end
end
end

%% 전이 함수 정의
function y = transition_function(x)
y = double(x >= 0 & x <= 1);
end

2.2 간소화된 EWT 구현 (초보자용)

function [modes, edges] = simplified_ewt(data, modeCount)
N = length(data);
dataFreq = fft(data);
mag = abs(dataFreq(1:floor(N/2)));

% 균등 분할 기반 경계 설정
binCount = floor(N/2);
stepSize = floor(binCount / modeCount);
edges = 1:stepSize:binCount;
edges(end) = binCount;

% 사각형 윈도우 기반 필터 생성
filters = zeros(modeCount, N);
for m = 1:modeCount
    startBin = edges(m);
    endBin = min(edges(m+1), N/2);
    filters(m, startBin:endBin) = 1;
    filters(m, N-endBin+2:N-startBin+2) = 1;
end

% 성분 추출
modes = zeros(modeCount, N);
for m = 1:modeCount
    resultFreq = dataFreq .* filters(m,:)';
    modes(m, :) = real(ifft(resultFreq));
end
end

3. 실용 예제

3.1 비정상 신호 분해

%% 다성분 신호 분석 예제
clear; clc; close all;

% 샘플링 설정
Fs = 1000;
t = 0:1/Fs:1.999;
T = length(t);

% 신호 구성 요소 생성
comp1 = 2.0 * sin(2*pi*5*t);                    % 저주파 성분
comp2 = 1.8 * sin(2*pi*(20 + 30*t).*t);          % 주파수 변조 성분
comp3 = 0.7 * sin(2*pi*110*t) .* exp(-((t-1)/0.15).^2); % 임펄스 성분
noise = 0.25 * randn(size(t));
compositeSignal = comp1 + comp2 + comp3 + noise;

% EWT 분석 수행
[imfs, limits] = empirical_wavelet_transform(compositeSignal, ...
    'numPeaks', 4, 'boundaryFactor', 0.35, 'waveletWidth', 0.09);

% 결과 시각화
figure('Position', [100, 100, 1300, 800], 'Color', 'white');

% 원본 신호
subplot(3, 4, 1);
plot(t, compositeSignal, 'Color', [0.2, 0.2, 0.2], 'LineWidth', 1.2);
xlabel('시간 (s)'); ylabel('진폭');
title('원본 신호'); grid on;

% 진폭 스펙트럼
subplot(3, 4, 2);
fAxis = (0:T-1)*(Fs/T);
plot(fAxis(1:T/2), abs(fft(compositeSignal)(1:T/2)), 'b', 'LineWidth', 1.2);
xlabel('주파수 (Hz)'); ylabel('진폭');
title('주파수 스펙트럼'); xlim([0, 150]); grid on;

% IMF별 추출 결과
labels = {'IMF1 (5Hz)', 'IMF2 (~35Hz)', 'IMF3 (110Hz)', '잔여 성분'};
colors = lines(4);

for k = 1:min(4, size(imfs,1))
    subplot(3, 4, 4+k);
    plot(t, imfs(k,:), 'Color', colors(k,:), 'LineWidth', 1.3);
    xlabel('시간 (s)'); ylabel('진폭');
    title(labels{k}); grid on;
end

% 재구성 검증
subplot(3, 4, 9);
recon = sum(imfs, 1);
plot(t, recon, 'k-', 'LineWidth', 1.4); hold on;
plot(t, compositeSignal, 'r--', 'LineWidth', 1);
legend('재구성', '원본', 'Location', 'best');
title('재구성 정확도'); grid on;

% 오차 분석
subplot(3, 4, 10);
errorSig = compositeSignal - recon;
plot(t, errorSig, 'Color', [0.5, 0, 0], 'LineWidth', 1.2);
xlabel('시간 (s)'); ylabel('오차');
title('재구성 오차'); grid on; ylim([-0.1, 0.1]);

% IMF별 주파수 분포
subplot(3, 4, 11);
for k = 1:size(imfs,1)
    spec = abs(fft(imfs(k,:)));
    plot(fAxis(1:T/2), spec(1:T/2), 'LineWidth', 1.1); hold on;
end
xlabel('주파수 (Hz)'); ylabel('진폭');
title('성분별 주파수 분포'); xlim([0, 150]); grid on;
legend('IMF1','IMF2','IMF3','IMF4');

% 에너지 기여도
subplot(3, 4, 12);
energies = sum(imfs.^2, 2);
pie(energies, sprintfc('IMF%d: %.1f%%', (1:length(energies))', energies'/sum(energies)*100));
title('에너지 기여도');

sgtitle('EWT 기반 신호 분해 결과', 'FontSize', 14, 'FontWeight', 'bold');

3.2 고장 진단 응용: 베어링 결함 탐지

%% 산업 설비 진단 예제
clear; clc; close all;

% 시뮬레이션 파라미터
Fs = 10000;
t = 0:1/Fs:0.6-1/Fs;

% 고장 특성 주파수 설정
innerRaceBPFI = 158;
ballPassFreq = 102;
resonanceFreq = 2200;

% 고장 관련 임펄스 생성
impulseTrain = zeros(size(t));
intervalInner = round(Fs / innerRaceBPFI);
for idx = 1:intervalInner:length(t)
    impulseTrain(idx) = 1.0;
end

% 공진 반응 모델링
decayRate = 600;
vibResponse = zeros(size(t));
for pos = 1:length(t)
    if impulseTrain(pos) == 1
        tau = t(pos:end) - t(pos);
        vibResponse(pos:end) = vibResponse(pos:end) + ...
            exp(-decayRate*tau) .* sin(2*pi*resonanceFreq*tau);
    end
end

% 배경 잡음 및 간섭 추가
backgroundNoise = 0.4 * randn(size(t));
powerLineInterference = 0.25 * sin(2*pi*60*t);
finalSignal = vibResponse + backgroundNoise + powerLineInterference;

% EWT 기반 분석
[imfs, bnds] = empirical_wavelet_transform(finalSignal, 'numPeaks', 5, 'boundaryFactor', 0.4);

% 고장 특징 추출
envelopeEnergy = zeros(size(imfs,1),1);
detectedFrequency = zeros(size(imfs,1),1);

for i = 1:size(imfs,1)
    analyticImf = hilbert(imfs(i,:));
    envelope = abs(analyticImf);
    envSpectrum = abs(fft(envelope));
    freqVec = (0:length(envSpectrum)-1)*Fs/length(envSpectrum);
    
    [~, peakIdx] = max(envSpectrum(1:round(300*length(envSpectrum)/Fs)));
    detectedFrequency(i) = freqVec(peakIdx);
    envelopeEnergy(i) = max(envSpectrum);
end

[~, dominantIMF] = max(envelopeEnergy);

% 시각화
figure('Position', [80, 80, 1500, 900], 'Color', 'w');

% 원본 진동 신호
subplot(3,5,1);
plot(t, finalSignal, 'k', 'LineWidth', 0.8);
xlabel('시간 (s)'); ylabel('가속도');
title('베어링 진동 신호'); grid on;

% 국부적 패턴 확대
subplot(3,5,2);
idxRegion = t>0.2 & t<0.25;
plot(t(idxRegion), finalSignal(idxRegion), 'r', 'LineWidth', 1.3);
title('임펄스 특성 확대'); grid on;

% 전체 주파수 분석
subplot(3,5,3);
fullSpec = abs(fft(finalSignal));
freqFull = (0:length(fullSpec)-1)*Fs/length(fullSpec);
plot(freqFull(1:length(fullSpec)/2), fullSpec(1:length(fullSpec)/2), 'b');
xlabel('주파수 (Hz)'); ylabel('진폭');
title('전체 스펙트럼'); xlim([0, 3000]); grid on;

% 포락선 스펙트럼 (원본)
subplot(3,5,4);
analyticOrig = hilbert(finalSignal);
envOrig = abs(analyticOrig);
specEnv = abs(fft(envOrig));
plot(freqFull(1:length(specEnv)/2), specEnv(1:length(specEnv)/2), 'g');
hold on; yline(innerRaceBPFI, 'r--', 'Label', '내륜 고장 주파수');
xlabel('주파수 (Hz)'); ylabel('진폭');
title('포락선 스펙트럼'); xlim([0, 300]); grid on;

% IMF 시각화 (상위 6개)
for k = 1:min(6,size(imfs,1))
    subplot(3,5,k+5);
    plot(t, imfs(k,:), 'LineWidth', 1.0);
    title(sprintf('IMF%d',k)); grid on;
end

% 주요 고장 성분의 포락선 분석
subplot(3,5,12);
selectedImf = imfs(dominantIMF,:);
analyticSel = hilbert(selectedImf);
envSel = abs(analyticSel);
specSel = abs(fft(envSel));
plot(freqFull(1:length(specSel)/2), specSel(1:length(specSel)/2), 'r', 'LineWidth', 1.5);
hold on; yline(innerRaceBPFI, 'r--', 'LineWidth', 1.5);
title(sprintf('IMF%d 포락선 스펙트럼', dominantIMF));
xlabel('주파수 (Hz)'); ylabel('진폭'); xlim([0, 300]); grid on;

% 진단 리포트
subplot(3,5,[14,15]);
axis off;
reportText = {
    '=== 지능형 진단 리포트 ===';
    '';
    sprintf('감지된 고장 주파수: %.1f Hz', detectedFrequency(dominantIMF));
    sprintf('결함 유형: %s', abs(detectedFrequency(dominantIMF) - innerRaceBPFI) < 5 ? ...
        '내륜 결함' : '볼 결함');
    '';
    sprintf('심각도 평가: %s', envelopeEnergy(dominantIMF) > 0.8 ? '높음' : '중간');
    '';
    '조치 권고:';
    '- 내륜 점검 필요';
    '- 조기 교체 권장';
    '- 모니터링 주기 단축';
};

text(0.1, 0.95, reportText, 'FontSize', 11, 'VerticalAlignment', 'top', ...
    'FontWeight', envelopeEnergy(dominantIMF) > 0.8 ? 'bold' : 'normal');
title('진단 결론', 'FontSize', 13, 'FontWeight', 'bold');

sgtitle('EWT 기반 베어링 고장 진단 시스템', 'FontSize', 15);

4. 성능 비교 및 평가

%% 다양한 기법의 성능 비교
function performance_comparison()
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
testSignal = sin(2*pi*10*t) + 0.7*sin(2*pi*45*t) + 0.5*sin(2*pi*110*t) + 0.15*randn(size(t));

% EWT 실행 시간 측정
tic;
[imfsEWT, ~] = empirical_wavelet_transform(testSignal, 'numPeaks', 3);
timeEWT = toc;

% EMD 실행 시간 측정
tic;
imfsEMD = emd(testSignal);
timeEMD = toc;

% 비교 지표 계산
orthogonalityEWT = calculate_cross_correlation(imfsEWT);
reconstructionErrorEWT = norm(testSignal - sum(imfsEWT,1)) / norm(testSignal);
spectralOverlapEWT = measure_spectral_interference(imfsEWT);

% 결과 시각화 생략 (실제 구현 시 추가)

fprintf('EWT 처리 시간: %.4f초\n', timeEWT);
fprintf('EMD 처리 시간: %.4f초\n', timeEMD);
end

function cc = calculate_cross_correlation(dataMatrix)
R = corrcoef(dataMatrix);
cc = sum(sum(triu(R,1))) / (size(dataMatrix,1)*(size(dataMatrix,1)-1)/2);
end

function overlap = measure_spectral_interference(dataMatrix)
overlap = 0;
for i = 1:size(dataMatrix,1)
    spec_i = abs(fft(dataMatrix(i,:)));
    for j = i+1:size(dataMatrix,1)
        spec_j = abs(fft(dataMatrix(j,:)));
        overlap = overlap + sum(min(spec_i, spec_j)) / (sum(spec_i) + 1e-10);
    end
end
end

5. 실용적 활용 가이드

5.1 파라미터 튜닝 전략

5.2 일반적인 문제 해결

• 성분 간 중첩 발생 시: boundaryFactor 감소 또는 numPeaks 증가
• 재구성 오차 큼: 필터 경계 확인, waveletWidth 조정
• 주요 성분 누락: 로그 스케일 적용 또는 피크 탐지 민감도 증가
• 처리 속도 느림: 신호 길이 단축 또는 간소화 버전 사용