What is wrong with my One Third Octave Filter for Frequency Weighted Acceleration

I am trying to obtain frequency weigted RMS acceleration through Matlab Code of One Third Octave Filter as per Method shown in the research work attached. While running the Matlab code it is showing the error, which I am not able to understand. What error my matlab code is having? The Matlab Code, Reference paper for code and Input Data (Acc. vs Time ) is attached.

The error obtained is:

??? Error using ==> butter at 41

Wn must be a one or two element vector.

Error in ==> One_Third_Octave_Filter at 59

[b,a]=butter(n,[wl wu]);

The Code is written here also:

%% One Third Octave Filter %%

clc;

clear all;

close all;

data = xlsread('Passive.xlsx');

Acc = data(:,2); % Passive sprung mass Acceleration vector

x_time= data(:,1); %Time vector

dt=0.001;

fs=1/dt; %Sampling Frequency%

% wc is vector of center frequencies of 1/3 octave bands, Hz %

wc=[1 1.25 1.6 2 2.5 3.15 4 5 6.3 8 10 12.5 16 20 25 31.5 40 50 63 80];

% weighting factor for centre frequency in 1/3 octave bands %

wi=[0.482 0.484 0.494 0.531 0.631 0.804 0.967 1.039 1.054 1.036 0.988 0.902 0.768 0.636 0.513 0.405 0.314 0.246 0.186 0.132];

%Lower Frequency band

Fl=(2^(-1/6))*wc;

%Higher Frequency band

Fh=(2^(1/6))*wc;

%Normalized Refernce Bandwidth

Br=((2^(1/6))-(2^(-1/6)));

% Reference Bandwidth quotient Qr

Qr=(1/Br);

%Designed Bandwidth quotiend Qd

n=4; %4th Order butter filter%

Qd=(((pi/(2*n))/sin(pi/(2*n)))*Qr);

%Set the ratio of the upper and lower band-edge frequencies to the mid-band frequency as ? or Beta,and it can be expressed by Qd %

Beta= ((1+sqrt(1+((2*Qd)^2)))/(2*Qd));

%The normalized cutoff frequencies are shown as follows:

wl=(wc/((fs/2)*Beta)); %fs is sampling frequency%

wu=((wc*Beta)/((fs/2)));

[b,a]=butter(n,[wl wu]);

y=filter(b,a,Acc);

Answer 1

[KALYAN ACHARJYA]

?

In the initial check, "wl" and "wu" are vector, it may be scalar data. Since it represents the lower and upper cuffs of frequencies in the Butterworth filter

>> whos wl
  Name      Size            Bytes  Class     Attributes
  wl        1x20              160  double              
>> whos wu
  Name      Size            Bytes  Class     Attributes
  wu        1x20              160  double   

Or

You may design multiple filters by passing individual single data (wl(i) and (wu)) one by one.

[Pankaj Sharma]

Dear KALYAN ACHARJYA

Shall I make this change to my code?

for i=1:length(wc);

[b(i),a(i)]=butter(n,[wl (i) wu(i)]);

Y(i)=filter(b,a,Acc);

end

will it resolve my problem?

[KALYAN ACHARJYA]

Sir, if you pass scalar (single value data) in cutoff frequencies, you can avoid coding error. Although I have no idea what exactly you want to do. Nevertheless, I did not examine any of the literature you attached.

for i=1:length(wu)
[b,a]=butter(n,[wl(i),wu(i)]);
freqz(b,a);
hold on;
end

It creats the multiple filters (20, as per length of wu or wl) , are you looking for that?

Please modify according to, so far you are clear why you got that error in your original code.

[Pankaj Sharma]

Thanks a lot for your valuable reply. Actually I am trying to prepare the code for Frequency Weighted RMS Acceleration in 1/3 Octave Frequency Band for the analysis of ride comfort as per ISO-2631-1. The paper I have attached is used as reference.

The output of this code through filter should be one single value (Scalar) know as RMS value, which otherwise is not coming as output and is coming out as vector.

[Pankaj Sharma]

Dear KALYAN ACHARJYA

I am facing another problem with another code. The code of a filter is prepared by Transfer Function as per the image attached (Image 1). The Image obtained by running Matlab code is a bit different when I see the X-axis value in terms of frequency, when compared with output image provided in Image 1. If you can suggest where I am making the mistake in my code.

Actually I am new to matlab so not very good in it.

The code use is as follows:

%% Weighting of Vertical acceleration as per ISO-2631-1 %%

num=[87.72 1138 11336 5453 5509];

den=[1 92.6854 2549.83 25969 81057 79783];

P=tf(num,den);

[h,f]= freqs(num,den);

F=(f/(2*pi));

mag=abs(h);

plot(F,mag);

xlabel('Frequency (Hz)','fontsize',14)

ylabel('Weighting Factor (Wk)','fontsize',14)

xlim([0 50])

Image-1 (Above)

Image Obtaiend by running above mentioned Matlab Code:

Image-2