Q Factor is not same

HI i have a matlab code but i am unable to find the mistake. i have attched the code.

Basically what i am unable to achieve is based on the Quality factor formula Q=f/f2-f1 we supposed to have the same quality factor in the input values alpha = 2*pi/Q. so, When, I change Amplitude or the frequency the Q factor goes far away from the right Q factor. In, unfortunately, I think I got some errors in the dB scale or something else. I must have Q factor in the output as the same Q factor in the input. If wouldn't mind I submit the Matlab code take look at and fix the problem. I cant figure out where is the error in the code that is messing up with this.

• Input:
• Digital Signal:
• F= 100 Hz
• Q= 500
• alpha= (f*pi)/Q
• A= 10

Result:

Q= fmax/f2-f1

Q= 237.5527

• Digital Signal:
• F= 100 Hz
• Q= 500
• A= 100

Result:

Q= 99.4541

Any help is appreciated. Thanks in advance.

clear all
clc
close all
fs= 10000
 
t = (0:1/fs:(1-(1/fs)));
f= 4000
alpha= (f*pi)/500
 y = 10*exp(-alpha*t).*sin(2*pi*f*t);
 figure
 plot(t,y)
 
 
 
N1=10000
S =fft(y,N1);
 
F1=S(1:N1/2); %half of spectra 
PF1=2*F1.*conj(F1)/(fs*N1); %Power spectrum density per 1 kHz
E=sum(PF1)%energy of signal V^2*T
Xdb = 20*log10(S);
LPF1=10*log10(PF1); %Power spectrum in dB scale
w=(1:N1/2);
LP=LPF1(1:N1/2);
w1=fs*w/N1;
figure 
plot(w1,LP)
 
 
BB= LP
CC= w1
E1 = max(LP)
t1= max(t)
talph1= E / -t1
 
%alpha1= max(BB*20*log(e))
 
[max_Z, max_index]=max(BB) %maximum value of impedance  
threedb=max_Z*sqrt(2)/2; %the 3db point 
[db_Z,index_db] = min(abs(BB.'-threedb)-max_Z)
                                            
%db3idx = zci(pwr-hpp);
[pk,loc] = findpeaks(BB,'Npeaks',1,'SortStr','descend');
%db3c= threedb-max_Z
ofst= 12
for k1 = 1:length(pk) 
 
    varmtx=[(BB(loc(k1)-ofst:loc(k1)));CC(loc(k1)-ofst:loc(k1));(BB(loc(k1):loc(k1)+ofst));CC(loc(k1):loc(k1)+ofst)]  
    dBpts(k1,1:2) = interp1(varmtx(1,:), varmtx(2,:), -(db_Z), 'linear','extrap');
    dBpts(k1,3:4) = interp1(varmtx(3,:),varmtx(4,:), -(db_Z), 'linear','extrap');
   
 
end 
 
figure(3)
plot(CC, BB)
hold on
for k1 = 1:length(pk)
    plot(dBpts(k1,:), -(db_Z), 'r+')
end
hold off
grid
 
F_1 = dBpts(k1,1:2)
F_1 = F_1(1)
F_2 = dBpts(k1,3:4)
F_2 =F_2(1)
BW  = (F_2)-(F_1) 
Q   = abs((loc-1)*(1/BW))

Answer 1

[David Goodmanson]

Hello Abdullah

be aware that what you are calling alpha is often called alpha/2, but I don't think that had anything to do with the problem you are having. I stayed with your definition of alpha.

Everything appears to be fine up until the width calculation. I don't know whether you intended that w denote circular frequency, but in this case it is (regular) frequency. To avoid ambiguity I went with frequency and did not use circular frequency. The code below finds the widths and the resulting Q.

clc
close all

fs= 10000;
t = (0:1/fs:(1-(1/fs)));
f= 4000;
alpha= (f*pi)/500;
 y = 10*exp(-alpha*t).*sin(2*pi*f*t);
 figure(1)
 plot(t,y); grid on
Q0 = pi*f/alpha
 
N1=10000;
S =fft(y,N1);
 
F1=S(1:N1/2); %half of spectra 
PF1=2*F1.*conj(F1)/(fs*N1); %Power spectrum density per 1 kHz
% E=sum(PF1);%energy of signal V^2*T
% Xdb = 20*log10(S);
LPF1=10*log10(PF1); %Power spectrum in dB scale
LP=LPF1(1:N1/2);

fvec=fs*(1:N1/2)/N1;   % frequency array
figure(2) 
plot(fvec,LP); grid on
xlim([3990 4010])

% find half power points by interpolation of LP
[LPmax ind0] = max(LP);
LPhalf = LPmax - 10*log10(2);            % half power points
indx = find(LP>LPhalf);
ind1 = [indx(1)-5:ind0];    % use a few more index points, 5 is somewhat arbitrary
f1 = interp1(LP(ind1),fvec(ind1),LPhalf,'spline');
ind2 = [ind0:indx(end)+5]; 
f2 = interp1(LP(ind2),fvec(ind2),LPhalf,'spline');
Q = f/(f2-f1)


  Q0 = 500.0000
  Q =  499.9917

[Abdullah Javed]

Thanks David