matlab, butterworth

Question

I need to apply a Butterworth filter. In order to prove my cutoff frequency, I wish to plot the residuals vs cutoff frequency, as shown in the image on the right:

Residuals are computed using the expression shown in the image on the left.

I want to apply a second-order low-pass Butterworth filter. My sample frequency is 40 Hz.

Using help butter reads:

"The cutoff frequency Wn must be 0.0 < Wn < 1.0, with 1.0 corresponding to half the sample rate".

This is my code, where x is the unfiltered data:

x = [386.321780100000;383.587759000000;381.800054600000;380.023872700000;376.516637700000;372.766842600000;367.938632100000;362.564671700000;356.119399400000;349.780736400000;345.535196800000;342.487775500000;343.327007400000;347.044262000000;352.155104000000;357.638520100000;363.112515900000;367.005736600000;369.969518400000;370.788505900000;370.727950800000;369.648907400000;371.655303500000;373.985894500000;374.662538300000;375.416536900000;377.716554500000;381.078738000000;384.060026700000;386.001162100000;386.230440100000;387.427096800000;388.211446400000;387.377444300000;387.621577500000;388.724233300000;389.431947300000;390.055687800000;389.897747100000;390.568256500000;390.976876100000;391.697054400000;392.307437700000;392.379471700000;393.675472100000;394.860379300000;395.811096400000;396.381388400000;398.085646700000;398.901758100000;400.678034200000;401.383476200000;402.159903200000;403.830729400000;404.584496200000;403.345648000000;402.534291300000;402.043024900000;400.591963700000;400.661529800000;401.214493200000;398.942282600000;396.136849800000;395.044892200000;391.983771800000;386.594017600000;383.401605600000;380.970551400000;376.310714500000;370.553521300000;366.878466900000;359.667106200000;352.458992900000;348.570689900000;343.701101100000;342.403147200000;345.523482600000;348.358285100000;351.430608400000;355.299118900000;358.338253800000;361.871167300000;365.253629000000;366.719133600000;371.786174300000;371.057137200000;371.100155400000;373.306458800000;376.272386700000;378.730711200000;382.095951400000;382.902277600000;383.461876500000;385.863290900000;386.271992900000;387.578084600000;388.651095700000;389.696072800000;390.282441000000;391.093495300000;391.289613800000;392.301607600000;392.159602900000;393.842911300000;393.546366500000;393.143725400000;394.750351600000;396.162576400000;397.561458700000;400.072731200000;401.171161400000;402.670562400000;403.090483400000;403.375118800000;403.742372900000;402.928992600000;403.625379400000;402.703021800000;404.445457900000;402.694859900000;404.492961300000;403.197389000000;400.246386500000;399.305927400000;395.243039500000;390.789534800000;386.215227800000;382.125630300000;378.740613300000;374.144220800000;369.667304900000;363.266125400000;358.880981700000;353.666361700000;352.930119300000;352.415712500000;354.932272100000;357.965910000000;363.210217400000;365.696508000000;368.583631400000;369.891869100000;370.457067900000;370.084533700000;371.254514800000;369.990608500000;372.512902900000;372.107388100000;372.389762900000;375.236476100000;377.180620200000;380.163700300000;383.752406800000;384.716569800000;386.690217000000;386.478825700000;387.136547300000;387.061047000000;387.614793500000;387.825040300000;388.859724400000;389.817438600000;390.783390500000;390.843715700000;391.629490100000;391.243866000000;391.117297600000;390.968419100000;392.245959800000;393.020590700000;394.630712100000;397.104264900000;397.023543100000;398.694405500000;399.038632200000;399.893340600000;401.723041400000;401.471017100000;402.798789400000;401.796557800000;401.510791300000;401.239454200000;402.555012900000;403.052252000000;402.218771700000;400.672248600000;398.567749700000;395.930037700000;393.687816300000;389.643818200000;385.300569800000;381.968732100000;379.096336600000;373.408813900000;368.723420200000;362.338487600000;357.422170600000;353.912378500000;351.497230900000;352.752737900000;355.645486800000;360.977540300000;365.115011900000;369.417573800000;374.421265000000;375.741261200000;379.482460500000;379.373812300000;381.356734700000;384.563742700000;385.998075700000;387.353239700000;387.846571200000;390.153373800000;394.202352800000;395.963324800000;398.282343000000;400.269374600000;400.557049800000;401.282810600000;402.133647600000;402.901177300000;404.506791700000;404.483256500000;405.463024500000;407.374795800000;407.037477400000;407.818079800000;408.022211200000;409.331968900000;410.195044600000;410.400538600000;411.138077700000;412.780554100000;413.371996900000;414.366633900000;416.607984600000;416.909471500000;417.903319900000;418.807274800000;418.131025900000;418.614937400000;419.355146400000;419.839171900000;420.622516800000;420.840936000000;421.415952700000;422.329106100000;419.791224700000;418.657891400000;416.303727400000;413.424186000000;409.239999100000;404.994917000000;401.405607300000;396.443609500000;391.830650800000;386.712834300000;379.305649600000;371.974335100000;367.876575200000;364.538574800000;364.350451800000;364.972675300000;368.720855400000;372.201352100000;374.288222600000;377.680482600000;378.799850200000;380.521041100000;380.363991200000;382.176218700000;384.164502100000;384.020557500000;387.386685800000;389.172806500000;391.352743100000;394.966380800000;397.758792800000;399.853885000000;399.654978000000;400.448029000000;402.805451900000;403.892553200000;403.862181700000;404.869756300000;406.689487700000;407.688509300000;407.899500500000;408.580087900000;410.215761300000;411.184311500000;412.070551500000;412.954576000000;414.074107800000;415.002936200000;416.180255400000;417.042582900000;417.210742100000;418.784222700000;419.417155900000;418.244356100000;418.892681200000;419.545868000000;419.460151500000;419.047418400000;419.924293500000;420.700464900000;420.588551400000;421.507250500000;422.020167300000;419.435987000000;417.067352800000;415.007963000000;411.769922000000;406.885412600000;402.982885100000;398.193117900000;390.465404700000;384.144538500000;375.178591500000;366.437035100000;358.028063500000;352.293817400000;349.424189400000;349.926651000000;352.430220700000;355.610311600000;360.417644500000;365.049574200000;369.087678900000;370.693150700000;371.469359300000;374.189527900000;376.566619700000;375.305055800000;375.891980400000;378.984007800000;382.437661200000;385.019999600000;387.613217500000;388.256001000000;390.440337200000;390.778502500000;394.132478900000;394.077478600000;395.934430300000;396.105431400000;396.078562100000;397.256080200000;397.439351600000;398.588204100000;398.132028000000;398.864585100000;400.041381600000;400.235601100000;401.558329500000;403.094163300000;403.243391900000;405.561711700000;405.504060100000;407.370466600000;408.555827900000;409.299491700000;410.153424800000;410.300822700000;410.367713700000;411.275381200000;410.871249900000;411.459785300000;412.533510600000;412.628470300000;414.186740100000;410.888678400000;411.371077500000;408.705744100000;406.903008300000;405.253466900000;400.188060400000;398.075866400000;394.070537500000;390.282135900000;386.964201200000;380.292136200000;374.835483700000;367.706128700000;365.908725400000;364.455562900000;364.563011500000;366.696664000000;370.685727600000;374.140337100000;375.529254600000;378.415772700000;379.025289000000;379.315148100000;379.507123000000;382.555575000000;385.180528900000;385.022058700000;388.587833500000;391.024301100000;395.008552400000;397.377134700000;399.993664600000;401.851884000000;403.349924700000;404.177882200000;404.530414400000;405.939468100000;405.167599800000;406.071639300000;406.301709800000;406.884170500000;408.251562800000;407.015528000000;407.044604300000;407.607805000000;408.133807700000;410.176208300000;410.295306300000;411.325515700000;412.468318200000;413.496999300000;414.096330200000;415.604314600000;416.924199100000;416.593331500000;418.376824200000;419.272729400000;418.891425800000;420.450437200000;419.941856000000;419.651679800000;421.103346200000;419.451771900000;417.418609800000;415.486484600000;413.071044300000;410.839635000000;405.672979600000;400.740956300000;395.880537300000;391.271580700000;385.603089900000;379.337494000000;373.292081100000;366.428944200000;359.433089600000;356.639741900000;356.870334800000;357.613582700000;360.618030500000;365.276770700000;367.955837500000;370.669455000000;374.535820600000;376.458963200000;377.488332300000;380.734952700000;380.967888300000;383.131782200000;382.965837100000;385.870830400000;388.718512600000;390.750714900000;394.412024500000;395.005689100000;396.315897800000;397.864541200000;399.388767400000;399.571308300000;400.029569200000;401.395641400000;402.391365000000;402.083481000000;402.528470400000;404.442816400000;404.884497900000;405.168726400000;405.465712800000;405.349603800000;406.533500000000;407.972460300000;409.293491400000;410.556036400000;412.133548900000;412.613902800000;412.779444300000;412.994332200000;413.556486900000;414.014261100000;415.313246400000;415.571968200000;416.646239500000;416.177362500000;416.654811500000;416.360201500000;414.128825500000;411.598404700000;410.069296900000;406.485724900000;402.166566900000;396.875429500000;393.467259500000;388.508421700000;381.249884300000;374.033839900000;367.407453700000;358.298466700000;353.014985400000;351.401607700000;351.559149200000;355.478413600000;361.579029700000;366.459723900000;370.482811700000;372.391993600000;374.133429600000;375.706849800000;376.255429100000;377.959570200000;378.150022100000;378.857670200000;381.393086700000;384.613589500000;387.134322900000;388.828588800000;388.884380500000;390.540899500000;391.078596600000;392.939800800000;393.131168700000;394.561911500000;394.940456400000;395.624279900000;396.631495200000;396.813280100000;397.271054100000;398.699067400000;397.960646900000;399.065647300000;400.495105800000;402.316133400000;402.820439400000;404.135953800000;404.902676700000;407.071295400000;408.901028400000;409.339043500000;410.856690000000;411.741902100000;411.572607400000;410.578750600000;410.866900900000;411.043940500000;412.270310500000;413.150341400000;413.186237200000;410.641084600000;407.405695900000;405.509177300000;402.693585500000;399.684704300000;394.426107300000;391.519533600000;388.224900900000;383.278789800000;378.307242300000;372.739602700000;368.478280700000;362.863358800000;361.304001000000;362.435018000000;364.634794200000;367.973272500000;370.867271600000;375.950784200000;377.760681400000;379.167124200000;379.646290600000;379.508894500000;379.045460000000;378.591352500000;384.507430200000;384.886844300000;383.932432300000;386.799684400000;391.060092500000;393.418588700000;394.132615800000;394.432453800000;398.324667600000;397.096586200000;398.005721400000;400.474732000000;400.328915300000;401.058421000000;403.238428000000;402.724730000000;403.253563800000;404.687953000000;405.079930400000;405.219330800000;407.047728400000;407.595529300000;408.057711800000;408.539252700000;409.836387000000;410.123266400000;410.313569400000;412.099500000000;412.793454600000;413.372957300000;415.066474200000;413.971619400000;414.254327400000;415.198047500000;416.134740700000;417.933600600000;417.211507700000;417.097779400000;417.708317200000;414.171904400000;412.326958700000;410.969281800000;403.879401900000;399.416247400000;393.297210900000;389.540637800000;384.152425200000;376.404144800000;372.813625800000;366.430311000000;361.093876700000;358.405677900000;358.336124000000;358.933857600000;363.550827300000;368.513013500000;373.555796400000;377.399377100000;382.541254000000;385.282657700000;389.342241000000;390.319845400000;390.834423800000;390.247456200000;392.237231500000;393.850930200000;398.146367200000;400.318802200000;404.104563600000;404.659164900000;];
wc = 0.001:0.001:1 % <----- !!!
R = zeros(1, length(wc)-1);
for nn = 1:length(wc)-1
    [B,A] = butter(2,wc(nn));
    xprime = filtfilt(B,A,x);
    R(1,nn) = sqrt(1/length(x)*(sum(abs(x-xprime).^2)))
end         
plot(wc(1:end-1)*20,R) % <------ !!!!

And I got this:

But this does not seem to be right to me. What am I doing wrong?

Answer 1

You need to set where you want your cutoff frequency. For example, if you know you want the -3dB point to be at 15 Hz, you can filter your data like so . .

fs = 40; %// 40 Hz as you said
fc = 15; %// your cutoff in Hz
[b,a] = butter(2, fc/(fs/2));

y = filter(b,a,x); %// output

Just change the value of fc to change the cutoff frequency . . .

Using filtfilt as you have used in your code will effectively double the order of your filter because you are filtering the signal twice (once forwards and once reversed). If you want zero-phase offline filtering with 2nd order, change the order to 1 before using filtfilt.

EDIT

To get the plot you are looking for given the information you have provided, I would use the following script ...

fs = 40;

x = [386.321780100000;383.587759000000;381.800054600000;380.023872700000;376.516637700000;372.766842600000;367.938632100000;362.564671700000;356.119399400000;349.780736400000;345.535196800000;342.487775500000;343.327007400000;347.044262000000;352.155104000000;357.638520100000;363.112515900000;367.005736600000;369.969518400000;370.788505900000;370.727950800000;369.648907400000;371.655303500000;373.985894500000;374.662538300000;375.416536900000;377.716554500000;381.078738000000;384.060026700000;386.001162100000;386.230440100000;387.427096800000;388.211446400000;387.377444300000;387.621577500000;388.724233300000;389.431947300000;390.055687800000;389.897747100000;390.568256500000;390.976876100000;391.697054400000;392.307437700000;392.379471700000;393.675472100000;394.860379300000;395.811096400000;396.381388400000;398.085646700000;398.901758100000;400.678034200000;401.383476200000;402.159903200000;403.830729400000;404.584496200000;403.345648000000;402.534291300000;402.043024900000;400.591963700000;400.661529800000;401.214493200000;398.942282600000;396.136849800000;395.044892200000;391.983771800000;386.594017600000;383.401605600000;380.970551400000;376.310714500000;370.553521300000;366.878466900000;359.667106200000;352.458992900000;348.570689900000;343.701101100000;342.403147200000;345.523482600000;348.358285100000;351.430608400000;355.299118900000;358.338253800000;361.871167300000;365.253629000000;366.719133600000;371.786174300000;371.057137200000;371.100155400000;373.306458800000;376.272386700000;378.730711200000;382.095951400000;382.902277600000;383.461876500000;385.863290900000;386.271992900000;387.578084600000;388.651095700000;389.696072800000;390.282441000000;391.093495300000;391.289613800000;392.301607600000;392.159602900000;393.842911300000;393.546366500000;393.143725400000;394.750351600000;396.162576400000;397.561458700000;400.072731200000;401.171161400000;402.670562400000;403.090483400000;403.375118800000;403.742372900000;402.928992600000;403.625379400000;402.703021800000;404.445457900000;402.694859900000;404.492961300000;403.197389000000;400.246386500000;399.305927400000;395.243039500000;390.789534800000;386.215227800000;382.125630300000;378.740613300000;374.144220800000;369.667304900000;363.266125400000;358.880981700000;353.666361700000;352.930119300000;352.415712500000;354.932272100000;357.965910000000;363.210217400000;365.696508000000;368.583631400000;369.891869100000;370.457067900000;370.084533700000;371.254514800000;369.990608500000;372.512902900000;372.107388100000;372.389762900000;375.236476100000;377.180620200000;380.163700300000;383.752406800000;384.716569800000;386.690217000000;386.478825700000;387.136547300000;387.061047000000;387.614793500000;387.825040300000;388.859724400000;389.817438600000;390.783390500000;390.843715700000;391.629490100000;391.243866000000;391.117297600000;390.968419100000;392.245959800000;393.020590700000;394.630712100000;397.104264900000;397.023543100000;398.694405500000;399.038632200000;399.893340600000;401.723041400000;401.471017100000;402.798789400000;401.796557800000;401.510791300000;401.239454200000;402.555012900000;403.052252000000;402.218771700000;400.672248600000;398.567749700000;395.930037700000;393.687816300000;389.643818200000;385.300569800000;381.968732100000;379.096336600000;373.408813900000;368.723420200000;362.338487600000;357.422170600000;353.912378500000;351.497230900000;352.752737900000;355.645486800000;360.977540300000;365.115011900000;369.417573800000;374.421265000000;375.741261200000;379.482460500000;379.373812300000;381.356734700000;384.563742700000;385.998075700000;387.353239700000;387.846571200000;390.153373800000;394.202352800000;395.963324800000;398.282343000000;400.269374600000;400.557049800000;401.282810600000;402.133647600000;402.901177300000;404.506791700000;404.483256500000;405.463024500000;407.374795800000;407.037477400000;407.818079800000;408.022211200000;409.331968900000;410.195044600000;410.400538600000;411.138077700000;412.780554100000;413.371996900000;414.366633900000;416.607984600000;416.909471500000;417.903319900000;418.807274800000;418.131025900000;418.614937400000;419.355146400000;419.839171900000;420.622516800000;420.840936000000;421.415952700000;422.329106100000;419.791224700000;418.657891400000;416.303727400000;413.424186000000;409.239999100000;404.994917000000;401.405607300000;396.443609500000;391.830650800000;386.712834300000;379.305649600000;371.974335100000;367.876575200000;364.538574800000;364.350451800000;364.972675300000;368.720855400000;372.201352100000;374.288222600000;377.680482600000;378.799850200000;380.521041100000;380.363991200000;382.176218700000;384.164502100000;384.020557500000;387.386685800000;389.172806500000;391.352743100000;394.966380800000;397.758792800000;399.853885000000;399.654978000000;400.448029000000;402.805451900000;403.892553200000;403.862181700000;404.869756300000;406.689487700000;407.688509300000;407.899500500000;408.580087900000;410.215761300000;411.184311500000;412.070551500000;412.954576000000;414.074107800000;415.002936200000;416.180255400000;417.042582900000;417.210742100000;418.784222700000;419.417155900000;418.244356100000;418.892681200000;419.545868000000;419.460151500000;419.047418400000;419.924293500000;420.700464900000;420.588551400000;421.507250500000;422.020167300000;419.435987000000;417.067352800000;415.007963000000;411.769922000000;406.885412600000;402.982885100000;398.193117900000;390.465404700000;384.144538500000;375.178591500000;366.437035100000;358.028063500000;352.293817400000;349.424189400000;349.926651000000;352.430220700000;355.610311600000;360.417644500000;365.049574200000;369.087678900000;370.693150700000;371.469359300000;374.189527900000;376.566619700000;375.305055800000;375.891980400000;378.984007800000;382.437661200000;385.019999600000;387.613217500000;388.256001000000;390.440337200000;390.778502500000;394.132478900000;394.077478600000;395.934430300000;396.105431400000;396.078562100000;397.256080200000;397.439351600000;398.588204100000;398.132028000000;398.864585100000;400.041381600000;400.235601100000;401.558329500000;403.094163300000;403.243391900000;405.561711700000;405.504060100000;407.370466600000;408.555827900000;409.299491700000;410.153424800000;410.300822700000;410.367713700000;411.275381200000;410.871249900000;411.459785300000;412.533510600000;412.628470300000;414.186740100000;410.888678400000;411.371077500000;408.705744100000;406.903008300000;405.253466900000;400.188060400000;398.075866400000;394.070537500000;390.282135900000;386.964201200000;380.292136200000;374.835483700000;367.706128700000;365.908725400000;364.455562900000;364.563011500000;366.696664000000;370.685727600000;374.140337100000;375.529254600000;378.415772700000;379.025289000000;379.315148100000;379.507123000000;382.555575000000;385.180528900000;385.022058700000;388.587833500000;391.024301100000;395.008552400000;397.377134700000;399.993664600000;401.851884000000;403.349924700000;404.177882200000;404.530414400000;405.939468100000;405.167599800000;406.071639300000;406.301709800000;406.884170500000;408.251562800000;407.015528000000;407.044604300000;407.607805000000;408.133807700000;410.176208300000;410.295306300000;411.325515700000;412.468318200000;413.496999300000;414.096330200000;415.604314600000;416.924199100000;416.593331500000;418.376824200000;419.272729400000;418.891425800000;420.450437200000;419.941856000000;419.651679800000;421.103346200000;419.451771900000;417.418609800000;415.486484600000;413.071044300000;410.839635000000;405.672979600000;400.740956300000;395.880537300000;391.271580700000;385.603089900000;379.337494000000;373.292081100000;366.428944200000;359.433089600000;356.639741900000;356.870334800000;357.613582700000;360.618030500000;365.276770700000;367.955837500000;370.669455000000;374.535820600000;376.458963200000;377.488332300000;380.734952700000;380.967888300000;383.131782200000;382.965837100000;385.870830400000;388.718512600000;390.750714900000;394.412024500000;395.005689100000;396.315897800000;397.864541200000;399.388767400000;399.571308300000;400.029569200000;401.395641400000;402.391365000000;402.083481000000;402.528470400000;404.442816400000;404.884497900000;405.168726400000;405.465712800000;405.349603800000;406.533500000000;407.972460300000;409.293491400000;410.556036400000;412.133548900000;412.613902800000;412.779444300000;412.994332200000;413.556486900000;414.014261100000;415.313246400000;415.571968200000;416.646239500000;416.177362500000;416.654811500000;416.360201500000;414.128825500000;411.598404700000;410.069296900000;406.485724900000;402.166566900000;396.875429500000;393.467259500000;388.508421700000;381.249884300000;374.033839900000;367.407453700000;358.298466700000;353.014985400000;351.401607700000;351.559149200000;355.478413600000;361.579029700000;366.459723900000;370.482811700000;372.391993600000;374.133429600000;375.706849800000;376.255429100000;377.959570200000;378.150022100000;378.857670200000;381.393086700000;384.613589500000;387.134322900000;388.828588800000;388.884380500000;390.540899500000;391.078596600000;392.939800800000;393.131168700000;394.561911500000;394.940456400000;395.624279900000;396.631495200000;396.813280100000;397.271054100000;398.699067400000;397.960646900000;399.065647300000;400.495105800000;402.316133400000;402.820439400000;404.135953800000;404.902676700000;407.071295400000;408.901028400000;409.339043500000;410.856690000000;411.741902100000;411.572607400000;410.578750600000;410.866900900000;411.043940500000;412.270310500000;413.150341400000;413.186237200000;410.641084600000;407.405695900000;405.509177300000;402.693585500000;399.684704300000;394.426107300000;391.519533600000;388.224900900000;383.278789800000;378.307242300000;372.739602700000;368.478280700000;362.863358800000;361.304001000000;362.435018000000;364.634794200000;367.973272500000;370.867271600000;375.950784200000;377.760681400000;379.167124200000;379.646290600000;379.508894500000;379.045460000000;378.591352500000;384.507430200000;384.886844300000;383.932432300000;386.799684400000;391.060092500000;393.418588700000;394.132615800000;394.432453800000;398.324667600000;397.096586200000;398.005721400000;400.474732000000;400.328915300000;401.058421000000;403.238428000000;402.724730000000;403.253563800000;404.687953000000;405.079930400000;405.219330800000;407.047728400000;407.595529300000;408.057711800000;408.539252700000;409.836387000000;410.123266400000;410.313569400000;412.099500000000;412.793454600000;413.372957300000;415.066474200000;413.971619400000;414.254327400000;415.198047500000;416.134740700000;417.933600600000;417.211507700000;417.097779400000;417.708317200000;414.171904400000;412.326958700000;410.969281800000;403.879401900000;399.416247400000;393.297210900000;389.540637800000;384.152425200000;376.404144800000;372.813625800000;366.430311000000;361.093876700000;358.405677900000;358.336124000000;358.933857600000;363.550827300000;368.513013500000;373.555796400000;377.399377100000;382.541254000000;385.282657700000;389.342241000000;390.319845400000;390.834423800000;390.247456200000;392.237231500000;393.850930200000;398.146367200000;400.318802200000;404.104563600000;404.659164900000;];
N = numel(x);

dw = 0.01;
wc = dw:dw:1-dw;

R=zeros(size(wc));
for nn =1:numel(wc)
    [b,a]=butter(1,wc(nn)); %//Potentially an order less because of filtfilt
    y = filtfilt(b,a,x);
    R(nn)=sqrt( sum((x-y).^2)/N );
end


fAxis = wc*(fs/2);

plot(fAxis,R);
xlabel('fc [Hz]'); ylabel('residuals')

If I use your data in my script, I get an output that looks like this. .

..which isn't too far of the output you posed in a link from your script ...