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

6 vues (au cours des 30 derniers jours)
Nathan Lively le 23 Juin 2022
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');
##### 12 commentairesAfficher 10 commentaires plus anciensMasquer 10 commentaires plus anciens
Mathieu NOE le 24 Juin 2022
ok but can you share some you already have and are meaningfull for you ?
Nathan Lively le 24 Juin 2022

Connectez-vous pour commenter.

### Réponse acceptée

Mathieu NOE le 28 Juin 2022
Nathan
as said above, here the full code :
%% playgroundFillGaps
clc;clear;close;format compact; format short g;
newFs = 96000; Nyq=newFs/2;
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 (°)');
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Plus de réponses (2)

Jon le 24 Juin 2022
Modifié(e) : Jon le 24 Juin 2022
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
##### 4 commentairesAfficher 2 commentaires plus anciensMasquer 2 commentaires plus anciens
Nathan Lively le 24 Juin 2022
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 le 24 Juin 2022
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.

Connectez-vous pour commenter.

Nathan Lively le 26 Juin 2022
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;
TF.Properties.VariableNames = {'frequency','magnitude','phase','coherence'};
%% Create timetable
TT = timetable('SampleRate',newFs);
TT.frequency = TF.frequency;
TT.IR = real(ifft(TT.Z,'symmetric')) * height(TT);
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));
%% Plot
semilogx(TF.frequency,TF.magnitude,TT.frequency,TT.magnitude2)
legend('original magnitude dB','fillsgaps magnitude dB','Location', 'Best','FontSize',30)
##### 9 commentairesAfficher 7 commentaires plus anciensMasquer 7 commentaires plus anciens
Nathan Lively le 27 Juin 2022
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 le 28 Juin 2022
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 !

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Digital Filtering dans Help Center et File Exchange

R2021a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by