0

高周波を除去するフィルターを実装したい。この例では、sin曲線を削除し、線形曲線を維持したいと思います。

編集：

コードを修正しましたが、フィルタリング用に実装した関数によってデータのエッジが大幅に変更され、受け入れられません。

clc; clear all;
xaxis = linspace(1, 10, 1000);
data = xaxis + sin(xaxis*3);
Nf = 2^12;
xAxisf = linspace(-10,10,Nf);
% plot(xaxis, data);

% FFT

xsize = numel(data);    
Xf = zeros([1 Nf]);
indices = Nf/2-floor(xsize/2):Nf/2-floor(xsize/2)+xsize - 1;
Xf(indices) = data;

% Xf = fftshift(Xf);
Xf = fft(Xf);
Xf = fftshift(Xf);

% plot 
Xfa = abs(Xf); plot(xAxisf, Xfa);

% generate super-gaussian filter function
Nf = numel(Xf);    
widthfilter = 0.12;
filterpower = 2;
filter = exp(-(xAxisf.^2./widthfilter^2).^filterpower);

% filter
filtertimes = 20;
Xf = Xf .* filter.^filtertimes;

% plot 
Xfa = abs(Xf); plot(Xfa);

% iFFt
Xfs = ifftshift(Xf);
Xif = ifft(Xfs);
% Xif = ifftshift(Xif);
result = abs(Xif);

plot(result(indices))
4

1 に答える 1

4

創刊：

    Xf = fftshift(data);     % NOT NEEDED
    Xf = fft(Xf);
    Xf = fftshift(Xf);

fft の前にデータを fftshift しないでください。シフトは fft の後にのみ必要です。これは、処理中に基数 n (おそらく 2) fft がデータを「デシメート」するためです。間引きされていないため、事前に修正する必要はありません。

2 番目の問題:

    Xfs = ifftshift(Xf);
    Xif = ifft2(Xfs);            
    Xif = ifftshift(Xif);   % NOT NEEDED

ifftshift は、入力として ifft が必要とするデータを再度デシメートします (fftshift を元に戻します)。入力が既に間引きされている場合、 ifft の出力は元の信号を再構成します。後にシフトシフトしないでください。

3 番目の問題:

    Xfs = ifftshift(Xf);
    Xif = ifft2(Xfs);       % USE ifft INSTEAD OF ifft2     
    Xif = ifftshift(Xiff);

なんでいきなり2D ifftに切り替えたの？

私はあなたのフィルターコードを詳しく見ていませんでしたが、ローパスフィルターが必要な場合は、中間点を中心に対称である必要があることに注意してください。そうしないと、周波数応答が対称的ではなく、多くの虚数になってしまいます。

そして、タイトルを変更してください。これは「フーリエ フィルター」ではありません。ウィンドウ法とfftを用いたローパスフィルターです。周波数空間でウィンドウを適用しているウィンドウ。

わかりました、遅くなり、前後に不機嫌になっています...ただ助けようとしています。私がコードを書いたほうが早いです。

コードでフィルターの効果を探している場合、フィルターのカットオフ周波数が高すぎるか、データ内の正弦波の周波数が低すぎるため、それができません。これは、入力正弦波の発振周波数を上げたバージョンです。

clc; clear all;
xaxis = linspace(1, 10, 1000);
data = xaxis + sin(xaxis*10);
% plot(xaxis, data);

% FFT
Xf = data;
Xf = fft(Xf);
Xf = fftshift(Xf);

% generate super-gaussian filter function
Nf = numel(data);
xAxisf = linspace(-5,5,Nf);
widthfilter = 0.1;
filterpower = 2;
filter = exp(-(xAxisf.^2./widthfilter^2).^filterpower);

% filter
filtertimes = 1;
Xf = Xf .* filter.^filtertimes;

% plot
Xfa = abs(Xf); plot(Xfa);

% iFFt
Xfs = ifftshift(Xf);
Xif = ifft(Xfs);
result = abs(Xif);

plot(result); hold on; plot(data,'r');
    legend('filtered','data');

就寝。おやすみ！私の公務を行いました:p

于 2013-01-25T10:30:40.420 に答える