MATLAB Answers

Q Factor is not same

24 views (last 30 days)
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
  1 Comment
Abdullah Javed
Abdullah Javed on 26 Feb 2021
Thanks David

Sign in to comment.

More Answers (0)

Tags

Products

Community Treasure Hunt

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

Start Hunting!