Problem with polyfit and polyval

Hello there, I took some manual measurements of the frequency response of a pre-amplifier and now I want to do a polyfit to get a smoother approximation (as I only measured some key points, the resolution is quite irregular) of the FR, which I do with polyfit. But so far, it returns me ZERO for all the polynomial coefficients! I'm probably doing a really simple mistake.

This is the code (the Amp.mat file is annexed):

%%Frequency Response Plot

close all; clear all

load('Amp.mat'); % Manual Measurements of the FR

Freq = Amp(:,1);          % Frequency
G = Amp(:,3)./Amp(:,2);   % Gain
Phase = Amp(:,4);         % Phase

GFit = polyfit(Freq,G,length(Freq));
GFitVal = polyval(GFit,Freq); 

PhaseFit = polyfit(Freq,Phase,length(Freq));
PhaseFitVal = polyval(PhaseFit,Freq); 

figure;

subplot(2,1,1)

semilogx(Freq,G, ...
         Freq,GFitVal); 
grid on;
legend('G', ...
       'GFitVal');
xlabel('Frequency (Hz)')
ylabel('Gain')
title('Gain x Frequency of the Pre-amplifier')

subplot(2,1,2)

semilogx(Freq,Phase, ...
         Freq,PhaseFitVal); 
grid on;
legend('Phase', ...
       'PhaseFit');
xlabel('Frequency (Hz)')
ylabel('Phase (degrees)')
title('Phase x Frequency of the Pre-amplifier')

% Impulse Response

x = 0:1e5;

GS = polyval(GFit,x);
PhaseS = polyval(PhaseFit,x);

FilterFreq = [flip(GS(2:end)).*exp(1i*(-flip(PhaseS(2:end)))) ...
              GS.*exp(1i*PhaseS)];

FilterIR = ifft(FilterFreq);

figure;

stem(FilterIR)
grid on;
xlabel('Samples (n)')
ylabel('Level')
title('Impulse Response of the Pre-amplifier')

Thanks in advance for the attention and possible help.

Answer 1

[Star Strider]

It took me a few minutes to come up with this purely signal processing approach, so no curve fitting required. Experiment with the numerator and denominator orders of the transfer function to get the result you want. (Another option would be the firls function, and there are other empirical methods to consider depending on what you want to do. There is also the entire System Identification Toolbox! You are doing system identification here.)

It should give you everything you want:

D = load('Philippe Amp.mat');
Freq = D.Amp(:,1);                                              % Hz Frequency Vector
Vi = D.Amp(:,3);
Vo = D.Amp(:,4);
H = Vo./Vi;                                                     % Amplitude Transfer Function
Phase = D.Amp(:,4);
W = Freq/max(Freq)*pi;                                          % Radian Frequency Vector

figure(1)                                                       % Look At The Data
subplot(2,1,1)
plot(Freq,Vi,  Freq,Vo)
grid
subplot(2,1,2)
plot(Freq, Phase)
grid

figure(2)
subplot(2,1,1)
plot(Freq, 20*log10(H))
ylabel('Magnitude (dB)')
grid
axis([xlim    -40  0])
subplot(2,1,2)
plot(Freq, Phase)
ylabel('Phase (degrees)')
grid

Hc = H .* exp(1i*pi*Phase/180);                                 % Create Complex Frequency Response
Hc(1) = interp1(Freq(2:end),Hc(2:end), 0, 'spline','extrap');   % Replace Initial ‘NaN’ Value

OrdNum = 3;
OrdDen = 5;
[b,a] = invfreqz(Hc, W, OrdNum, OrdDen);                        % Transfer Function Coefficients

figure(3)
freqz(b, a, 4096)                                               % Bode Plot

figure(4)
impz(b, a)                                                      % Impulse Response

Answer 2

[John D'Errico]

You start off doing a bad thing in these two lines:.

GFit = polyfit(Freq,G,length(Freq));
PhaseFit = polyfit(Freq,Phase,length(Freq));

Warning: Polynomial is not unique; degree >= number of data points. 
> In polyfit (line 70) 
Warning: Polynomial is not unique; degree >= number of data points. 
> In polyfit (line 70)

You cannot estimate more coefficients than you have data points. If you try, you will get garbage. You will NEVER get a "smoother" approximation using that approach.

Worse:

size(Freq)
ans =
    63     1

This is insane! Do you honestly think that you can estimate that high a degree a polynomial? Whoever told you to do so was flat out wrong. The result will be completely meaningless garbage. Oh yeah, that is what you got.

IMHO, polyfit should be limited to never more than a 5th degree polynomial or so. But even then, some problems will fail to allow you to estimate even a quadratic polynomial. Some times, you might get something vaguely intelligent with a degree as high as 10. So sadly, there is no constraint in polyfit, besides the common sense of the user, and their understanding of what they are doing. But remove common sense and any understanding, and you should expect meaningless results.

[Philippe]

Ok, my bad, so indeed I was doing something really wrong. Nobody told me to do that, it was something I wanted to test. Do you have any suggestions as to how to estimate an approximation of the measured FR function for me to get the IR with the ifft?

[John D'Errico]

I've often said this. Polynomials are nice. They are basically Taylor series, and we use series for just about everything, right? They teach us about series in early calc classes, so they must be good. And if a low order polynomial works, then a high order one must be better!

WRONG. High order polynomials are bad. Essentially, always. You will get garbage out from this computation.

Something I should have done before, and something you should ALWAYS learn to do before you build any models, is plot the data. Think about what you see. Apply common sense.

plot(Freq,G)

plot(Freq,Phase)

Do either of those curves look like anything a polynomial would fit? (Hint: NO!)

Why not? What shape do they have at zero? You appear to have virtually a singularity. Do polynomials represent functions with singularities well? NO.

So it matters not what order polynomial you tried here, you would expect garbage.

If you are willing to try to smooth through the initial apparent singularity, use a LOW order polynomial.

P3 = polyfit(Freq/1e6,Phase,3)
P3 =
     1.4667      -9.3682      -19.745       149.63

Note that here, I scaled Freq, dividing by 1e6. Otherwise, you will force the polynomial to cube numbers on the order of 1e6, and then perform arithmetic with them. It will generate garbage again.

plot(Freq,Phase,'r.')
hold on
plot(Freq,polyval(P3,Freq/1e6),'g-')

This points out that you don't really have a singularity at zero. Just a function that is RAPIDLY varying. Do polynomials represent that well? NO!

In fact, I'm surprised that polyfit produced non-garbage here. Go to a much higher order, and it would have started to produce garbage quickly though.

With a polynomial model, stay low order. If you need a high order polynomial, then you are doing something wrong.

[Philippe]

You're right again, as this is a FR of a filter, I believe a bandpass filter function like Butterworth or the Raised Cossine would be more adequate, do you know any function that would approximate my data points to these types of filters?

[John D'Errico]

A good idea when your function varies so rapidly at one end of the domain is to consider a transformation. A simple test of that fact is to try this:

semilogx(Freq,Phase,'r.')

That suggest a log transformation makes sense. So, model the responses as a function of the log of Freq.

As well, even the above curve will NOT be terribly amenable to a polynomial fit. While it would not be as obscene as your original fit might have been, you would still need to use a fairly high order polynomial to have some chance at a decent fit.

The point is, this is the time when a least squares spline makes sense. Using my SLM toolbox (on the file exchange), I get this curve in the domain of log(Freq).

slm = slmengine(log(Freq),Phase,'knots',12,'plot','on')

In the original domain, we would see this:

plot(Freq,Phase,'r.',Freq,slmeval(log(Freq),slm),'g-')

But it is a spline. You cannot look at the coefficients and make any sense of them, nor can you even write them down in any intelligent form.

[Image Analyst]

I'm still wondering why it needed smoothing in the first place. Looks pretty smooth to me even without filtering. Why is smoothing even needed???

Answer 3

[Image Analyst]

What order polynomial do you want to fit it to? When you do this:

GFit = polyfit(Freq,G,length(Freq));

you're trying to fit your data to a 63rd order polynomial, which is insane to begin with, but the error message says that you can't fit a 63rd order polynomial to an array with only 63 points. You'd need at least 64 points for that.

Perhaps you want to use a sliding polynomial of order 2 or 3 like you can with sgolayfit().