Q Factor is not same

12 views (last 30 days)
Abdullah Javed
Abdullah Javed on 25 Feb 2021
Commented: Abdullah Javed on 26 Feb 2021
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))

Accepted Answer

David Goodmanson
David Goodmanson on 26 Feb 2021
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

More Answers (0)

Tags

Products

Community Treasure Hunt

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

Start Hunting!