0

Matlab で単純なフィルターの逆フーリエ変換を見つけようとしています。最初のケース (sinc フィルター/「ブリック ウォール」) では、関数を使用してifft、t = 0 を中心とする sinc である時間領域関数を見つけます。

ここで、単純なチェビシェフ フィルターの時間領域関数を見つけたいと思います。しかし、何らかの理由で、以下のコードでは、時間軸が正しくないインパルス応答が得られるようです。t = 0 を中心とする同様の見た目の sinc 関数を期待すべきではありませんか?

fo = 10; %frequency of the sine wave
Fs = 100; %sampling rate
Ts = 1/Fs; %sampling time interval
t = -1+Ts:Ts:1-Ts; %sampling period
freq = -Fs/2:(Fs/length(t)):Fs/2;

%% Sinc with bandwidth = fo. This works!
y = 0.5*sinc(2*fo*t);
YfreqDomain1 = fft(y);
figure('Name','Brick wall sinc filter (freq)');
plot(freq(1:length(y)),2/length(y)*fftshift(abs(YfreqDomain1)))
y_ret1=ifft(YfreqDomain1,'nonsymmetric');
figure('Name','Brick wall sinc filter (time)');
plot(t,y_ret1);

%% Chebyshev with bandwidth fo. This gives me a strange result. 
[b,a] = cheby1(6,0.1,2*fo/Fs);  % 6th order, 0.1dB ripple
[YfreqDomain2 w] = freqz(b,a,length(t),'whole');
figure('Name','Chebyshev Filter (freq)');
plot(freq(1:length(YfreqDomain2)), 2/length(y)*fftshift(abs(YfreqDomain2)));
figure('Name','Chebyshev Filter (time)');
y_ret2=ifft(YfreqDomain2,'nonsymmetric');
plot(t,y_ret2);
4

1 に答える 1

2

計算には 2 つの問題があります。

最初に、非常に狭い時間ウィンドウで時間領域フィルター係数を評価します。これにより、フィルターが打ち切られ、その特性が変化します。

第 2 に、時間領域ベクトルと周波数領域ベクトルのどのインデックスがそれぞれ時間点 0 と周波数 0 に対応するかを適切に追跡していません。これは別の方法で行うことができます。ここでは、常に と を選択しt(1) = 0、プロットのみf(1) = 0に適用することにしました。fftshift

あなたのコードの私の修正版は次のとおりです。

fo = 10;        % frequency of the sine wave
Fs = 100;       % sampling rate
Ts = 1 / Fs;    % sampling time interval
n = 1000;       % number of samples

% prepare sampling time vector such that t(1) = 0
t = (0 : n - 1)';
t = t - n * (t >= 0.5 * n);
t = t / Fs;

% prepare frequency vector such that f(1) = 0;
f = (0 : n - 1)' / n;
f = f - (f >= 0.5);
f = f * Fs;


%% sinc filter with bandwidth fo

% define filter in time domain
s = 2*fo/Fs * sinc(2*fo*t);

% transform into frequency domain
Hs = fft(s);
% since the filter is symmetric in time, it is purely real in the frequency
% domain. remove numeric deviations from that:
Hs = real(Hs);

subplot(2, 4, 1)
plot(fftshift(t), fftshift(s))
ylim([-0.1 0.25])
title('sinc, time domain')

subplot(2, 4, 2)
plot(fftshift(f), real(fftshift(Hs)), ...
    fftshift(f), imag(fftshift(Hs)))
ylim([-1.1 1.1])
title({'sinc, frequency domain', 'real / imaginary'})

subplot(2, 4, 3)
plot(fftshift(f), abs(fftshift(Hs)))
ylim([-0.1 1.1])
title({'sinc, frequency domain', 'modulus'})


%% Chebyshev filter with bandwidth fo

% define filter in frequency domain
[b,a] = cheby1(6, 0.1, 2*fo/Fs);  % 6th order, 0.1 dB ripple
Hc = freqz(b, a, n, 'whole', Fs);

% transform into time domain
c = ifft(Hc);

% determine phase such that phase(1) is as close to 0 as possible
phase = fftshift(unwrap(angle(fftshift(Hc))));
phase = phase - 2*pi * round(phase(1) /2/pi);

subplot(2, 4, 5)
plot(fftshift(t), fftshift(c))
title('Chebyshev, time domain')
ylim([-0.1 0.25])

subplot(2, 4, 6)
plot(fftshift(f), real(fftshift(Hc)), ...
    fftshift(f), imag(fftshift(Hc)))
ylim([-1.1 1.1])
title({'Chebyshev, frequency domain', 'real / imaginary'})

subplot(2, 4, 7)
plot(fftshift(f), abs(fftshift(Hc)))
ylim([-0.1 1.1])
title({'Chebyshev, frequency domain', 'modulus'})

subplot(2, 4, 8)
plot(fftshift(f), fftshift(phase))
title({'Chebyshev, frequency domain', 'phase'})

結果は次のとおりです。

ご覧のとおり、sinc フィルターとチェビシェフ フィルターは、周波数応答のモジュラスに関しては似ていますが、位相に関しては大きく異なります。その理由は、チェビシェフ フィルターが因果フィルターであるためです。つまり、時間領域の係数は t < 0 の場合は 0 に制限されます。これは、実際の物理システムに実装されるフィルターの自然な特性です。

于 2014-12-09T20:36:36.423 に答える