How to fill missing frequency response data at the beginning and end?

fillmissing and interp1 work great for gaps in the middle of an array, but I can't find anything that works well for filling the missing data at the beginning and end. It seems like maybe the Curve Fitting toolbox might do it, but I have searched around and can't find any examples of someone filling missing data at the beginning and end.

PROBLEM: Some of the frequency response data I get from an audio analyzer does not include all frequency bins from 0Hz to Nyquist.

SOLUTION: The best solution I have found is to estimate the magnitude at 0Hz and then use fillmissing for the data in between. The results are not great. I could easily make a better estimation by just draing with my eyes.

Here's an example where fillmissing was used to fill the data between 0-10Hz and 20-48kHz. I have tried all of the other methods (spline, etc.) and the all deliver much worse results. Thanks for any suggestions!!

S_TF = fillmissing(S_TF,'linear','EndValues','extrap');

[Mathieu NOE]

hello

seems to me there is no "scientific" possible approach to that problem

it just a matter of what data / plot you want to see and implement the code that realise it

simply decide what should be the amplitude at beginning and end points and trace a line , polynomial , exponential of what you prefer....

[Nathan Lively]

Thanks Mathieu. Well, I know there is a solution because I see it being done is some other audio analyzers. There's an app called Crosslite, for example. If I import the data there, it's fills in the beginning and end quite nicely. I want to know how it is accomplished.

[Mathieu NOE]

hello Nathan

found that : crosslite Archives - Sound Design Live with you as the author...so you are in audio equalization stuff - that's funny because that was also a topic for me a while ago even though I never pratice it at home (jounger I was a bit like an audio fanatic but it has calm down with age)

Ok - I don't know crosslite works in the details and how they fill the missing ends , but if by chance you had just a simple picture of one case treated by crosslite maybe we can figure out how to implement it in matlab.

[Nathan Lively]

Yes, that's me!

Well, I just ran a test and the results are less surprising. I guess Crosslite is not doing anything special. Still wish I could do something better.

[Mathieu NOE]

why not simply make a linear extrapolation based on low end / high end slope ?

[Adam Danz]

I agree with @Mathieu NOE's initial statement. fillmissing uses neighboring points to interpolate values and your end points do not have sufficient neighboring points. Futhermore, you should not rely on interpolating polynomials to extrapolate values at the end points.

https://ece.uwaterloo.ca/~dwharder/NumericalAnalysis/06LeastSquares/extrapolation/complete.html

If you have an expected function your data should follow, then you might have a shot and extracting meaningful exterpolation of the data.

[Nathan Lively]

Thanks @Mathieu NOE

I'd love to do that, but I don't know how to consistently tell it where to calculate the slope. In this example above, it would be from about 10-80Hz, but in others it will be slightly different. It will always be where the response starts to slope down from right to left, which is easy to see with my eyes, but I don't know hot to tell MATLAB to look for a downward trend.

[Nathan Lively]

Thanks @Adam Danz!

[Mathieu NOE]

@Nathan Lively

you have some data to share ?

[Nathan Lively]

@Mathieu NOE I have...a lot of data. :)

I actually have an online database the exchange of loudspeaker reference data called Tracebook. You are required to create a login, but after that, it's free for anyone to use. Under every graph there are links to download the data.

[Mathieu NOE]

ok but can you share some you already have and are meaningfull for you ?

[Nathan Lively]

@Mathieu NOE sure!

EAW KF300e. Complete data.

https://trace-book.org/measurement/eaw-kf300e-full-range-passive

https://s3.amazonaws.com/appforest_uf/f1642268846106x688929829017490800/KF300e%20on-axis.csv

RCF TTL 4-A. Complete data.

https://trace-book.org/measurement/rcf-ttl-4-a-c1

https://s3.amazonaws.com/appforest_uf/f1628873958066x749579241040526000/RCF.TTL4A.C1.2.74m%20-%20Cmplx.7150wCC.csv

NEXO PS10R2. Incomplete data

https://trace-book.org/measurement/nexo-ps10r2-35-ps10r2-wb

https://s3.amazonaws.com/appforest_uf/f1638022388563x970089448815184100/PS10R2%20WB.csv

Answer 1

[Mathieu NOE]

Nathan

as said above, here the full code :

%% playgroundFillGaps
clc;clear;close;format compact; format short g;
newFs = 96000; Nyq=newFs/2;
% Load data
TF = readmatrix('K3-110-0.99m-NR2.txt',"NumHeaderLines",9); % {'frequency','magnitude','phase','coherence'};
freq = TF(:,1);
mag_dB = TF(:,2);
phas_deg = TF(:,3);
coherence = TF(:,1);
%keep only data below Fs/2 freq
ind = (freq<=newFs/2);
freq = freq(ind);
mag_dB = mag_dB(ind);
phas_deg = phas_deg(ind);
coherence = coherence(ind);
% add dummy DC value (replicate first available data)
freq = [0; freq];
mag_dB = [mag_dB(1); mag_dB];
phas_deg = [phas_deg(1); phas_deg];
coherence = [coherence(1); coherence];
% make both ends have natural behaviour  = 
% high pass filter like behaviour in low freq end (below 20 Hz)
% low pass filter like behaviour in high freq end (above 20 kHz)
% apply frequency window = butterworth band pass filter
f_low = 20;      % lower cut off freq Hz
f_high = 20e3;      % higher cut off freq Hz
N_bpf = 2;           % filter order
[b,a] = butter(N_bpf,2/newFs*[f_low f_high]);
frf_filter = freqz(b,a,freq,newFs);
frf = 10.^(mag_dB/20).*exp(j*pi/180*phas_deg); % initial FRF
frf_filtered = frf.*frf_filter;% complex initial FRF + bandpass window FRF
%% Plot
figure(1)
subplot(2,1,1),semilogx(freq,mag_dB,freq,20*log10(abs(frf_filtered)))
subplot(2,1,2),semilogx(freq,phas_deg,freq,180/pi*angle(frf_filtered))
%% fft method
frf = frf_filtered;
% create fliiped negative freq FRF data (conjuguate)
if mod(length(frf),2)==0 % iseven
    frf_sym = conj(frf(end:-1:2));
else
    frf_sym = conj(frf(end-1:-1:2));
end

fir = real(ifft([frf; frf_sym])); % negative and positive frequency frf converted into IR
% fir = fir(1:10000); % truncation is possible if fir decays enough
frfid = freqz(fir,1,freq,newFs);
%% plot
figure(2)
subplot(311),plot(0,0,1:length(fir),fir);
legend('no data','identified FIR model (fft)');
xlabel('samples');
ylabel('amplitude');
subplot(312),semilogx(freq,20*log10(abs(frf)),freq,20*log10(abs(frfid)));
ylim([max(mag_dB)-80  max(mag_dB)+10]);
legend('input model','identified FIR model (fft)');
xlabel('frequency (Hz)');
ylabel('FRF modulus (dB)');
subplot(313),semilogx(freq,180/pi*angle(frf),freq,180/pi*angle(frfid));
legend('input model','identified FIR model (fft)');
xlabel('frequency (Hz)');
ylabel('FRF angle (°)');

Answer 2

[Jon]

You could fit your frequency response data with a transfer function model, and then use the model to generate the missing data. You might have to put in some constraints on your fitting to match what you "know" about the steady state (zero frequency) gain of the system is.

If you can make some assumptions (you have some knowledge) of the characteristics of the device at low frequency, you could just make a transfer function model for that with some matching condition at zero and the start of your data

Either way, as the others have pointed out, unless you have some physical basis for making the model, there is no correct way of generating the missing data, the behavior could be anything where you don't have data

[Nathan Lively]

Thanks @Jon.

> fit your frequency response data with a transfer function model

Does that mean: Compare data with a known good model where I have all of the data from beginning to end?

I don't know how to do that. Should I search for "MATLAB model comparison"?

[Jon]

If you have the System Identification Toolbox, this gives some background on how you would fit a model to frequency response data

https://www.mathworks.com/help/ident/gs/estimate-models-using-frequency-domain-data.html

[Nathan Lively]

Cool. Looks like misdata might be the way to go.

https://www.mathworks.com/help/ident/ug/handling-missing-data-and-outliers.html

[Nathan Lively]

Hey Jon, I don't know if you have any bandwidth to demo what you were talking about with some of the measurements I posted above. Otherwise, I'll check it out this weekend.

Answer 3

[Nathan Lively]

Ok, I had an idea to use fillgaps on the impulse response generated by ifft. I spent a few hours on it this weekend, but I can't seem to get it to work. I have tried several variations, but I think I've reached the limits of my understanding of FFT. Please let me know if you see any errors.

I kept reading about how hard it is to fill gaps at the beginning and end of arrays so I thought it would be a novel solution to convert the frequency domain data back to the time domain, circshift the impulse response so that the gap would be in the middle instead of the ends, then use fillgaps. Again, not sure why it didn't work.

%% playgroundFillGaps
clc;clear;close;format compact; format short g;
newFs = 96000; Nyq=newFs/2;

% Load data
url1 = 'https://www.sounddesignlive.com/wp-content/uploads/2022/06/K3-110-0.99m-NR2.txt';
TF = readtable(url1,'ReadVariableNames',false);
TF.Properties.VariableNames = {'frequency','magnitude','phase','coherence'};

%% Create timetable
TT = timetable('SampleRate',newFs);
TT.frequency = TF.frequency;
TT.Z = [db2mag(TF.magnitude) .* exp(1j*(deg2rad(TF.phase)))];
TT.IR = real(ifft(TT.Z,'symmetric')) * height(TT);

% Add a gap
TT.IRgap = TT.IR;
TT.IRgap(1:21) = NaN;

%% Fill gap
TT.IRshifted = circshift(TT.IRgap,Nyq);
TT.IRgapFill = fillgaps(TT.IRshifted);
TT.IRgapFill = circshift(TT.IRgapFill,Nyq*-1);

%% Back to frequency domain
TT.Z2 = fft(TT.IRgapFill) / height(TT);
TT.magnitude2 = mag2db(abs(TT.Z2));
TT.phase2 = rad2deg(angle(TT.Z2));

%% Plot
semilogx(TF.frequency,TF.magnitude,TT.frequency,TT.magnitude2)
legend('original magnitude dB','fillsgaps magnitude dB','Location', 'Best','FontSize',30)

[Mathieu NOE]

hello again

seems here data has been added below 10 Hz and above 20 kHz. This extra data does not really make sense of what a mesured loudpeaker would deliver

should be more like the red strips I plotted and next topic is how I will do it in matlab

mag plot above (in dB) , phase below.

[Nathan Lively]

hey @Mathieu NOE, yes, it's all noise down there, but my understanding is that it is necessary to have all of the data from 0Hz to Nyquist to accurately recreate the impulse response and move back and forth between the the time domain and frequency domain. That being said, yes, I'd much rather have higher SNR and have the loudpseaker response continue lower instead of running into the noise floor.

[Nathan Lively]

I realized this morning that I am probably doing this wrong because I don't quite undertand the relationship between the data in the frequency domain and where it appears in the time domain. I'll try to make some new code that make the gap in the frequency domain instead of the time domain.

[Nathan Lively]

Ok, this time I tried putting the gap in the frequency data. Results are still incorrect. Still unsure what the problem is.

%% playgroundFillGaps
clc;clear;close;format compact; format short g;
newFs = 96000; Nyq=newFs/2;

% Load data
url1 = 'https://www.sounddesignlive.com/wp-content/uploads/2022/06/K3-110-0.99m-NR2.txt';
TF = readtable(url1,'ReadVariableNames',false);
TF.Properties.VariableNames = {'frequency','magnitude','phase','coherence'};

% Add a gap
TF.magnitudeGap = TF.magnitude; TF.phaseGap = TF.phase;
TF.magnitudeGap(1:21) = NaN; TF.phaseGap(1:21) = NaN;

%% Create timetable
TT = timetable('SampleRate',newFs);
TT.frequency = TF.frequency;
TT.Z = [db2mag(TF.magnitude) .* exp(1j*(deg2rad(TF.phase)))];
TT.Zgap(:) = NaN; TT.Zgap(22:end) = [db2mag(TF.magnitudeGap(22:end)) .* exp(1j*(deg2rad(TF.phase(22:end))))];
TT.IR = real(ifft(TT.Z,'symmetric')) * height(TT);
TT.IRgap(:) = NaN; TT.IRgap(22:end) = real(ifft(TT.Z(22:end),'symmetric')) * height(TT.IRgap(22:end));

%% Fill gap
TT.IRshifted = circshift(TT.IRgap,Nyq);
TT.IRgapFill = fillgaps(TT.IRshifted);
TT.IRgapFill = circshift(TT.IRgapFill,Nyq*-1);

%% Back to frequency domain
TT.Z2 = fft(TT.IRgapFill) / height(TT);
TT.magnitude2 = mag2db(abs(TT.Z2));
TT.phase2 = rad2deg(angle(TT.Z2));

%% Plot
nexttile
semilogx(TF.frequency,TF.magnitude,TT.frequency,TT.magnitude2)
legend('original magnitude dB','fillsgaps magnitude dB','Location', 'Best','FontSize',30)
nexttile
plot(TT.Time,TT.IR,TT.Time,TT.IRgapFill)
xlim([seconds(0) seconds(0.01)])
legend('original ifft','ifft gap filled','Location', 'Best','FontSize',30)

[Mathieu NOE]

hello again

yes , if you want "perfect" IR , you need a transfer function measured from DC to Fs/2

if not the reconstructed IR by ifft will look distorded

[Nathan Lively]

I am getting closer, I think. This time I apply circshift to the known good data first, then fillgaps. I'm still missing something, though.

%% playgroundFillGaps
clc;clear;close;format compact; format short g;
newFs = 96000; Nyq=newFs/2;

% Load data
url1 = 'https://www.sounddesignlive.com/wp-content/uploads/2022/06/K3-110-0.99m-NR2.txt';
TF = readtable(url1,'ReadVariableNames',false);
TF.Properties.VariableNames = {'frequency','magnitude','phase','coherence'};

% Add a gap
TF.magnitudeGap = TF.magnitude; TF.phaseGap = TF.phase;
TF.magnitudeGap(1:21) = NaN; TF.phaseGap(1:21) = NaN;

%% Create timetable
TT = timetable('SampleRate',newFs);
TT.frequency = TF.frequency;
TT.Z = [db2mag(TF.magnitude) .* exp(1j*(deg2rad(TF.phase)))];
TT.Zgap(:) = NaN; TT.Zgap(22:end) = [db2mag(TF.magnitudeGap(22:end)) .* exp(1j*(deg2rad(TF.phase(22:end))))];
TT.IR = real(ifft(TT.Z,'symmetric')) * height(TT);
TT.IRgap(:) = NaN; TT.IRgap(22:end) = real(ifft(TT.Z(22:end),'symmetric')) * height(TT.IRgap(22:end));

%% Fill gap
TT.IRshifted(:) = TT.IRgap;
[~,irPeakIdx] = max(abs(TT.IRshifted(22:end))); % Find peak
shiftAmount = irPeakIdx * -1 +1 + round(height(TT.IRshifted(22:end))/2); % Calculate shift
TT.IRshifted(22:end) = circshift(TT.IRshifted(22:end),shiftAmount);
TT.IRgapFill = fillgaps(TT.IRshifted);
TT.IRgapFill = circshift(TT.IRgapFill,shiftAmount*-1 -21);

%% Back to frequency domain
TT.Z2 = fft(TT.IRgapFill) / height(TT);
TT.magnitude2 = mag2db(abs(TT.Z2));
TT.phase2 = rad2deg(angle(TT.Z2));

%% Plot
nexttile
semilogx(TF.frequency,TF.magnitude,TT.frequency,TT.magnitude2)
legend('original magnitude dB','fillsgaps magnitude dB','Location', 'Best','FontSize',30)
nexttile
plot(TT.Time,TT.IR,TT.Time,TT.IRgapFill)
xlim([seconds(0) seconds(0.001)])
legend('original ifft','ifft gap filled','Location', 'Best','FontSize',30)

[Mathieu NOE]

hello again

tried to play with my code below -

we can "tweak" the data to have 0 to Fs/2 completely filled

now the make sure both ends of the frequency range behaves like it shoud, I applied a butterworth filter as a frequency weighting window, to force the gain to be decreasing to zero as freq goes down until DC.

same idea for the higher freq limit

then I do the ifft stuff and we an see that the IR with all 96000 samples gives perfect match with the windowed FRF

but i's a very long IR, so I looked at what happens if you shorten it :

if we keep only the first 10,000 samples , there will be some approximations :

see the line :

% fir = fir(1:10000); % truncation is possible if fir decays enough

here the full code :

%% playgroundFillGaps
clc;clear;close;format compact; format short g;
newFs = 96000; Nyq=newFs/2;
% Load data
TF = readmatrix('K3-110-0.99m-NR2.txt',"NumHeaderLines",9); % {'frequency','magnitude','phase','coherence'};
freq = TF(:,1);
mag_dB = TF(:,2);
phas_deg = TF(:,3);
coherence = TF(:,1);

%keep only data below Fs/2 freq
ind = (freq<=newFs/2);
freq = freq(ind);
mag_dB = mag_dB(ind);
phas_deg = phas_deg(ind);
coherence = coherence(ind);

% add dummy DC value (replicate first available data)
freq = [0; freq];
mag_dB = [mag_dB(1); mag_dB];
phas_deg = [phas_deg(1); phas_deg];
coherence = [coherence(1); coherence];

% make both ends have natural behaviour  = 
% high pass filter like behaviour in low freq end (below 20 Hz)
% low pass filter like behaviour in high freq end (above 20 kHz)
% apply frequency window = butterworth band pass filter

f_low = 20;      % lower cut off freq Hz
f_high = 20e3;      % higher cut off freq Hz
 N_bpf = 2;           % filter order
[b,a] = butter(N_bpf,2/newFs*[f_low f_high]);
frf_filter = freqz(b,a,freq,newFs);

frf = 10.^(mag_dB/20).*exp(j*pi/180*phas_deg); % initial FRF

frf_filtered = frf.*frf_filter;% complex initial FRF + bandpass window FRF

%% Plot
figure(1)
subplot(2,1,1),semilogx(freq,mag_dB,freq,20*log10(abs(frf_filtered)))
subplot(2,1,2),semilogx(freq,phas_deg,freq,180/pi*angle(frf_filtered))

%% fft method
frf = frf_filtered;

% create fliiped negative freq FRF data (conjuguate)
if mod(length(frf),2)==0 % iseven
    frf_sym = conj(frf(end:-1:2));
else
    frf_sym = conj(frf(end-1:-1:2));
end
    
fir = real(ifft([frf; frf_sym])); % negative and positive frequency frf converted into IR
% fir = fir(1:10000); % truncation is possible if fir decays enough
frfid = freqz(fir,1,freq,newFs);

%% plot
figure(2)
subplot(311),plot(0,0,1:length(fir),fir);
legend('no data','identified FIR model (fft)');
xlabel('samples');
ylabel('amplitude');
subplot(312),semilogx(freq,20*log10(abs(frf)),freq,20*log10(abs(frfid)));
ylim([max(mag_dB)-80  max(mag_dB)+10]);
legend('input model','identified FIR model (fft)');
xlabel('frequency (Hz)');
ylabel('FRF modulus (dB)');
subplot(313),semilogx(freq,180/pi*angle(frf),freq,180/pi*angle(frfid));
legend('input model','identified FIR model (fft)');
xlabel('frequency (Hz)');
ylabel('FRF angle (°)');

[Nathan Lively]

Thanks @Mathieu NOE, this great!

Interesting, this is very similar to another function I wrote called noiseReduction. I guess I just thought it would be better to fill in any missing data with as much accuracy as possible before doing the filtering, but as you have pointed out, it's just noise down there anyway. Thanks again!

I'm not sure how to accept your answer above, since it is a comment.

[Mathieu NOE]

hello Nathan

glad my work culd be of some use !

I copied it as an aswer if you want to accepted , that'd great !

all the best for the future !