MATLAB: Fourier Graphs of the lowpass and high pass filtered signals are looking weird

fourier transformsignal processing

Hello everyone,

I am currently doing a project in MATLAB and I am using a lowpass filter in a signal with a cutoff frequency of 30Hz and then I find the fourier transform of it and I graph it. I want to know if the actual graph is right or am I doing something wrong? I am going to attach the part of the code for the filters and plots.

This is the plot that I get for the low pass filter:

It shouldn't go to 0 after 30 Hz? or around 30 Hz?

And for the high pass filter:

I get it to be 0 around until like 190 Hz, but why is it not 270Hz or closer? I did use smaller frequencies for the HPF and they were more accurate until when they stayed 0.

%LPF 30Hz
y_LPF_30 = lowpass(y,30,Fs);
yf_LPF_30 = fft(y_LPF_30);
% 30Hz
%set up

%% Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L.
P2_LPF30 = abs(yf_LPF_30/L);
P1_LPF30 = P2_LPF30(1:L/2+1);
P1_LPF30(2:end-1) = 2*P1_LPF30(2:end-1);
figure (5)
%plot

subplot(3,2,1);
plot(f,P1_LPF30) 
title('Single-Sided Amplitude Spectrum of LPF 30Hz')
xlabel('f (Hz)')
ylabel('|P1(f)|')
%HPF 270Hz
y_HPF_270 = highpass(y,270,Fs);
yf_HPF_270 = fft(y_HPF_270);
%270Hz
%set up
P2_HPF270 = abs(yf_HPF_270/L);
P1_HPF270 = P2_HPF270(1:L/2+1);
P1_HPF270(2:end-1) = 2*P1_HPF270(2:end-1);
%plot
subplot(3,2,6);
plot(f,P1_HPF270) 
title('Single-Sided Amplitude Spectrum of HPF 270Hz')
xlabel('f (Hz)')
ylabel('|P1(f)|')

Best Answer

Last update:

Running your Lowpass code:

%LPF 30Hz
y_LPF_30 = lowpass(y,30,Fs);
yf_LPF_30 = fft(y_LPF_30);
% 30Hz
%set up
%% Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L.
P2_LPF30 = abs(yf_LPF_30/L);
P1_LPF30 = P2_LPF30(1:L/2+1);
P1_LPF30(2:end-1) = 2*P1_LPF30(2:end-1);
figure
%plot
plot(f,P1_LPF30) 
title('Single-Sided Amplitude Spectrum of LPF 30Hz')
xlabel('f (Hz)')
ylabel('|P1(f)|')

gives this output:

Update:

For the low-pass try adding some properties such as:

y_LPF_30 = lowpass(y,30,Fs,'StopbandAttenuation',100);

Also, please show us the plots generated by this:

lowpass(y,30,Fs,'StopbandAttenuation',100);

lowpass(___) with no output arguments plots the input signal and overlays the filtered signal.

Hope it helps.

Earlier:

The filters look ok. some example data would be helpful to debug.

Also, it might just be a typo, but you wrote: P2_HPF270 = abs(yf_HPF_220/L);

Related Solutions

MATLAB: How to make the high pass filter more accurate

If my cutoff frequency is 270 Hz, why do I get still frequencies below 270?

Unless it's an ideal highpass filter, most filter implementations show a transition band where the attenuation starts to take effect. You can find more details here:

https://www.mathworks.com/help/signal/ref/highpass.html#mw_7e643f9d-3d99-43c8-afb5-585aede4cdcd

To get steeper response, try this:

highpass(y,270,Fs,'Steepness',0.99,'StopbandAttenuation',100)

MATLAB: 50 Hz interference and noise reduction from ECG.

hello

If you want only to remove the 50 Hz tone , use a notch filter as shown in example below. I used your data and guessed the sampling frequency from the fact that the peak was at 50 Hz.

So sampling at 10 kHz is pretty high for an ECG signal, you could easily go down by factor 10.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
signal = importdata('Sinal1.txt');
dt = 1e-4; % 0.1 milli seconds
samples = length(signal);
Fs = 1/dt; % sampling frequency (Hz)
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = Fs;    % 
Overlap = 0.75;
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
%% notch filter section %%%%%%
% H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
fc = 50;       % notch freq
wc = 2*pi*fc;
Q = 7;         % adjust Q factor for wider (low Q) / narrower (high Q) notch
                % at f = fc the filter has gain = 0
w0 = 2*pi*fc/Fs;
alpha = sin(w0)/(2*Q);
            b0 =   1;
            b1 =  -2*cos(w0);
            b2 =   1;
            a0 =   1 + alpha;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha;
            
% analog notch (for info)
num1=[1/wc^2 0 1];
den1=[1/wc^2 1/(wc*Q) 1];
% digital notch (for info)
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
freq = linspace(fc-1,fc+1,200);
[g1,p1] = bode(num1,den1,2*pi*freq);
g1db = 20*log10(g1);
[gd1,pd1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
gd1db = 20*log10(gd1);
figure(1);
plot(freq,g1db,'b',freq,gd1db,'+r');
title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');
legend('analog','digital');
xlabel('Fréquence (Hz)');
ylabel(' dB')
 
% now let's filter the signal 
signal_filtered = filtfilt(num1z,den1z,signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display : averaged FFT spectrum before / after notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq,fft_spectrum] = myfft_peak(signal, Fs, NFFT, Overlap);
sensor_spectrum_dB = 20*log10(fft_spectrum);% convert to dB scale (ref = 1)

[freq,fft_spectrum_filtered] = myfft_peak(signal_filtered, Fs, NFFT, Overlap);
sensor_spectrum_filtered_dB = 20*log10(fft_spectrum_filtered);% convert to dB scale (ref = 1)
figure(2),semilogx(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid on
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before notch filter','after notch filter');
xlabel('Frequency (Hz)');ylabel(' dB')
function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples)) = signal;
    signal = s_tmp;
    samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows
    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select);
freq_vector = (select - 1)*Fs/nfft;
end

Best Answer

Related Solutions

MATLAB: How to make the high pass filter more accurate

MATLAB: 50 Hz interference and noise reduction from ECG.

Related Question