matlab - 1/3 octave spectra from fft

I have got a signal in frequency domain. This is a frequency response function from software, so I can do nothing about it and have to leave it in frequency domain. I want to transfer the data to 1/3 octave spectra and I need some help on this.

In my opinion, I just need a couple of band pass filters which in the frequency domain is simply by using the starting and ending limit of each band to cut the data. In other words, I ignore all other data points out of a band region.

To get the band spectra, originally I just calculate the square of each fft point inside a band and sum them before calculating the square root. However, the result is different with direct filtering in the time domain, which I tried using a random signal. It makes sense to me because the 1/3 octave spectra should be independent of the total fft points. Obviously by simply calculating the sum of squares of fft components, the value of power in a band increases with the increase of the points included.

Then I tried to explore the basic theory underlining the fft based filters. This is where I got very confused. I do not know which part of FT frequency response function belongs to? It is not DFT in my opinion and not DFS either because the whole process of transferring motion equation is based on continuous FT except the frequency response function is calculated at a series of discontinuous points. My understanding is that each component of 1/3 octave spectra is a square root of the total power in that band. Therefore possibly I can get something from Parseval's theorem. However based on continuous FT transform, I could not figure out how to derive the power in a frequency range.

I would like to know the way to implement the filter and underlying theory if possible. Please give me some suggestions or alternatively some references that I can refer to.

Thanks a lot.