低通滤波器实现正确还是错误?
Low pass filter implementation Correct or wrong?
我一直在研究 2 个测量旋转轴振动的传感器信号。由于信号中存在残余噪声。我试图通过去除趋势、零填充和应用低通滤波器来过滤它。下面我附上了滤波前后的信号图。过滤后的信号变化很大,这让我思考我是否真的以正确的方式进行操作。
我的 Matlab 代码是
X = xlsread(filename,'F:F');
Y = xlsread(filename,'G:G');
%Calculate frequency axis
fs = 1e6 ; % Sampling frequency (Hz)
NFFT = 2^nextpow2(length(X)); % Zero padding to nearest N power 2
df = fs/NFFT;
dt = 1/df;
%Frequency Axis defintion
f = (-(fs-df)/2:df:(fs-df)/2)';
X(2^ceil(log2(length(X))))=0;
Y(2^ceil(log2(length(Y))))=0;
%calculate time axis
T = (dt:dt:(length(X)*dt))';
subplot(2,2,1)
plot(T,X);
xlabel('Time(s)')
ylabel('X amplitude')
title('X signal before filtering')
subplot(2,2,2)
plot(T,Y);
xlabel('Time(s)')
ylabel('Y amplitude')
title('Y signal before filtering')
X = detrend(X,0); % Removing DC Offset
Y = detrend(Y,0); % Removing DC Offset
% Filter parameters:
M = length(X); % signal length
L = M; % filter length
fc = 2*(38000/60); % cutoff frequency
% Design the filter using the window method:
hsupp = (-(L-1)/2:(L-1)/2);
hideal = (2*fc/fs)*sinc(2*fc*hsupp/fs);
h = hamming(L)' .* hideal; % h is our filter
% Zero pad the signal and impulse response:
X(2^ceil(log2(M)))=0;
xzp = X;
hzp = [ h zeros(1,NFFT-L) ];
% Transform the signal and the filter:
X = fft(xzp);
H = fft(hzp)';
X = X .* H;
X = ifft(X);
relrmserrX = norm(imag(X))/norm(X); % checked... this for zero
X = real(X)';
% Zero pad the signal and impulse response:
Y(2^ceil(log2(M)))=0;
xzp = Y;
hzp = [ h zeros(1,NFFT-L) ];
% Transform the signal and the filter:
Y = fft(xzp);
H = fft(hzp)';
Y = Y .* H;
Y = ifft(Y);
relrmserrY = norm(imag(Y))/norm(Y); % check... should be zero
Y = real(Y)';
我画了过滤后的图,你可以在图片中看到有明显的偏差。我只想过滤噪声,但信号似乎会丢失其他成分,如果那是正确的过滤方法,我会感到有点困惑。
任何建议、提示或想法都会有所帮助。
最后我想绘制 X 与 Y 的关系图以给出轴振动的轨道。另请在下面找到另一张未过滤和过滤轨道的图片。正如您在图片中看到的那样,原始轨道也发生了变化(左图有很多噪音)。
P.S.: 我没有DSP工具箱
你的FFT和IFFT没有问题
您可以使用频率非常低的简单正弦波来测试您的代码。最终输出应该是(几乎)相同的正弦波,因为你正在对它进行低通滤波。
您已将 X(和 Y)定义为从 0 到某个值。但是,您已将 H 从 - (L-1)/2 定义为某个正值。从数学上讲,这很好,但是您只是简单地采用 H 的 fft。当您将 ffts 相乘时,Matlab 认为这与 X 具有相同的时间尺度!
因此,实际上,您已经采用了 XHfft(delta(t-d)) 的 fft,其中 d 是产生的时移。您可以在频域或时域中撤消此时移。
我一直在研究 2 个测量旋转轴振动的传感器信号。由于信号中存在残余噪声。我试图通过去除趋势、零填充和应用低通滤波器来过滤它。下面我附上了滤波前后的信号图。过滤后的信号变化很大,这让我思考我是否真的以正确的方式进行操作。
我的 Matlab 代码是
X = xlsread(filename,'F:F');
Y = xlsread(filename,'G:G');
%Calculate frequency axis
fs = 1e6 ; % Sampling frequency (Hz)
NFFT = 2^nextpow2(length(X)); % Zero padding to nearest N power 2
df = fs/NFFT;
dt = 1/df;
%Frequency Axis defintion
f = (-(fs-df)/2:df:(fs-df)/2)';
X(2^ceil(log2(length(X))))=0;
Y(2^ceil(log2(length(Y))))=0;
%calculate time axis
T = (dt:dt:(length(X)*dt))';
subplot(2,2,1)
plot(T,X);
xlabel('Time(s)')
ylabel('X amplitude')
title('X signal before filtering')
subplot(2,2,2)
plot(T,Y);
xlabel('Time(s)')
ylabel('Y amplitude')
title('Y signal before filtering')
X = detrend(X,0); % Removing DC Offset
Y = detrend(Y,0); % Removing DC Offset
% Filter parameters:
M = length(X); % signal length
L = M; % filter length
fc = 2*(38000/60); % cutoff frequency
% Design the filter using the window method:
hsupp = (-(L-1)/2:(L-1)/2);
hideal = (2*fc/fs)*sinc(2*fc*hsupp/fs);
h = hamming(L)' .* hideal; % h is our filter
% Zero pad the signal and impulse response:
X(2^ceil(log2(M)))=0;
xzp = X;
hzp = [ h zeros(1,NFFT-L) ];
% Transform the signal and the filter:
X = fft(xzp);
H = fft(hzp)';
X = X .* H;
X = ifft(X);
relrmserrX = norm(imag(X))/norm(X); % checked... this for zero
X = real(X)';
% Zero pad the signal and impulse response:
Y(2^ceil(log2(M)))=0;
xzp = Y;
hzp = [ h zeros(1,NFFT-L) ];
% Transform the signal and the filter:
Y = fft(xzp);
H = fft(hzp)';
Y = Y .* H;
Y = ifft(Y);
relrmserrY = norm(imag(Y))/norm(Y); % check... should be zero
Y = real(Y)';
我画了过滤后的图,你可以在图片中看到有明显的偏差。我只想过滤噪声,但信号似乎会丢失其他成分,如果那是正确的过滤方法,我会感到有点困惑。 任何建议、提示或想法都会有所帮助。
最后我想绘制 X 与 Y 的关系图以给出轴振动的轨道。另请在下面找到另一张未过滤和过滤轨道的图片。正如您在图片中看到的那样,原始轨道也发生了变化(左图有很多噪音)。
P.S.: 我没有DSP工具箱
你的FFT和IFFT没有问题
您可以使用频率非常低的简单正弦波来测试您的代码。最终输出应该是(几乎)相同的正弦波,因为你正在对它进行低通滤波。
您已将 X(和 Y)定义为从 0 到某个值。但是,您已将 H 从 - (L-1)/2 定义为某个正值。从数学上讲,这很好,但是您只是简单地采用 H 的 fft。当您将 ffts 相乘时,Matlab 认为这与 X 具有相同的时间尺度!
因此,实际上,您已经采用了 XHfft(delta(t-d)) 的 fft,其中 d 是产生的时移。您可以在频域或时域中撤消此时移。