Effacer les filtres
Effacer les filtres

Problem in filtering my desired frequency

10 vues (au cours des 30 derniers jours)
Abdullah Alasfour
Abdullah Alasfour le 20 Juil 2023
Commenté : Paul le 29 Juil 2023
Hello guys, I have this problem with filtering, that is supposed to allow the frequency 730 to pass while block the other frequencies 440 and 490. The problem it didn't pass my required frequency 730
Based on the instruction given to us, I was supposed to build I deal filter that is constant around my desired frequency X(w) = 1. So, I did inverse fourier transform ( h(t) ) with paper, and it is integral from 2*720*pi to 2*740*pi and got this " ht= abs( (1/2*pi*j*t).*( exp(j*1480*pi*t) - exp(j*1440*pi*t))); " and plotted it in matlab. Any help :') ?
dt= 0.0001;
t = 0:dt:10;
xt = 4*cos(2*pi*730*t) +3*cos(2*pi*440*t) + 5*cos(2*pi*490*t);
subplot(2,2,1)
plot(t,xt)
axis([ 0 0.02 -20 20 ])
title(' Ploting x(t)')
xlabel('time')
ylabel('x(t)')
grid on
ht= abs( (1/2*pi*j*t).*( exp(j*1480*pi*t) - exp(j*1440*pi*t)));
subplot(2,2, 2)
plot(t,ht)
title(' plotting h(t)')
axis( [0 20 0 70] )
xlabel(' Time')
ylabel(' h(t)')
grid on
yt= conv(xt , ht, 'same');
subplot(2,2, 3)
plot(t,yt)
title(' plotting y(t)')
axis( [0 20 0 100] )
xlabel(' Time')
ylabel(' y(t)')
grid on
fs = 1/dt; %sampling rate sample/sec
N = length(t); % number of samples
f = ( 0: N-1)*(fs/N);
Yw= abs(fft( yt ));
subplot(2,2,4)
plot( f , Yw)
title(' Ploting Y(w)')
xlabel('frequency')
ylabel('Y(w)')
axis([ 0 3000 0 1.5*max(Yw)])
grid on

Réponse acceptée

Paul
Paul le 21 Juil 2023
Modifié(e) : Paul le 22 Juil 2023
Hi Abdullah,
A couple of issues to address in the anlaysis.
I think you started by trying to define an ideal filter in the frequency domain based on its Continuous Time Fourier Transform (CTFT), like so.
syms w t real
H(w) = rectangularPulse(2*sym(pi)*720,2*sym(pi)*740,w);
Its impulse response is
h(t) = rewrite(ifourier(H(w),w,t),'exp');
[num,den] = numden(h(t));
h(t) = num/den
h(t) = 
The code in the question defined h(t) as the magnitude (I'm not sure why) of this signal
( (1/2*sym(pi)*1j*t).*( exp(1j*1480*sym(pi)*t) - exp(1j*1440*sym(pi)*t)))
ans = 
Apparently there were some missing parentheses in the denominator of the leading term, added here
( (1/(2*sym(pi)*1j*t)).*( exp(1j*1480*sym(pi)*t) - exp(1j*1440*sym(pi)*t)))
ans = 
which matces h(t).
However, the ideal band pass filter has to also account for negative frequencies. So we take H(w) as defined above and add its mirror image.
H(w) = H(w) + H(-w);
Now the impulse response is
h(t) = ifourier(H(w),w,t);
[num,den] = numden(h(t));
h(t) = num/den
h(t) = 
h(t) is basically the sum of two sinc signals, a fact that will come into play below. Also, the impulse response is non-causal, i.e., it's non-zero for t < 0, and symmetric around t = 0, facts that will come into play below. Also, it has infinite duration, i.e., it never decays exactly to zero, which will also be important.
Define x(t) and its CTFT
x(t) = 4*cos(2*sym(pi)*730*t) + 3*cos(2*sym(pi)*440*t) + 5*cos(2*sym(pi)*490*t)
x(t) = 
X(w) = fourier(x(t),t,w)
X(w) = 
The CTFT of the output is then
Y(w) = X(w)*H(w);
And the output is
y(t) = simplify(ifourier(Y(w),w,t))
y(t) = 
As expected, the ideal bandpass filter exactly recovers the 730 Hz component and eliminates the other two.
If we want to compute the output of the filter via the convolution integral, it would look like this.
%syms tau real
%int(x(tau)*h(t-tau),tau,-inf,inf)
I commented it out because the computation was exceeding the Answers run time limit.
Now, if we want to try to do the analysis numerically in discrete time via the FFT, we have to deal with the fact that h(t), and x(t) for that matter, is infinite duration, and non-causal and symmetric, because FFT only applies to finite duration signals. The former means we'll have to pick a window over which h(t) is large enough to matter. The latter means we need to have our time vector be symmetric around t = 0, extending to negative time, as opposed to starting at t = 0 as in the Question. For x(t), we'll also have to pick a window as well.
Looking at h(t), it falls off as 1/(t*pi), so I'll make the h(t) window be -5 <= t <= 5, i.e., 1/(5*pi) seems pretty small. I'll arbitrarily use the same window for x(t).
Assuming the 10kHz sampling frequency, the sampled time vector is
Ts = 1e-4;
nx = -50000:50000;
nh = nx;
tvals = nx*Ts;
Evaluate x(t) and h(t) at the sampling times.
xf = matlabFunction(x(t));
hf = matlabFunction(h(t));
xval = xf(nx*Ts);
hval = hf(nh*Ts);
Recall that h(t) is really the sum of sinc functions, so we have to manually deal with the divide by 0 at t = 0 (or rewrite h(t) in terms of sinc I suppose).
limit(h(t),t,0)
ans = 
40
hval(nh==0) = 40;
Now we can evaluate samples of y(t) using the convolution sum. We need to multiply it by Ts to have it approximate the convolution integral, because h(t) and x(t) are both finite signals (i.e., don't include impulses).
yval = conv(hval,xval)*Ts; % uses 'full' as the default
If we wanted to plot the values of y in the time domain, the corresponding time vector would be
tyval = Ts*((0:(numel(yval)-1)) + nx(1) + nh(1));
tyval covers a larger range than the time window used for h(t) and x(t), which is a consequence of convolution.
tyval([1 end])
ans = 1×2
-10 10
In the frequency domain we take the fft of the input samples and multiply by Ts so the result approximates the CTFT of x(t)
X = fft(xval)*Ts; fx = (0:numel(X)-1)/numel(X)/Ts;
Do the same for the samples of the impulse response
H = fft(hval)*Ts; fh = fx;
Plot the frequency responses X and H, limited to the frequencies less than 1000 Hz for clarity
figure
plot(fx,abs(X),fh,abs(H)),xlim([0 1000]),xlabel('Hz');ylabel('Amplitude')
legend('Input','Filter')
The filter has nice, near rectangular pulse at 730 Hz. The input has the spikes at the three frequencies.
The amplitude of X at 730 Hz is ~20, whereas X(w) above has a Dirac delta at that frequency.
Recall that in the discrete-time analysis we're not working with x(t), but actually a windowed version of x(t). The CTFT of the 730 Hz component of the windowed signal is
Xr(w) = fourier(4*cos(2*sym(pi)*730*t)*rectangularPulse(-5,5,t),t,w)
Xr(w) = 
And its amplitude at 730 Hz is
limit(Xr(w),2*sym(pi)*730)
ans = 
20
as seen in the plot.
Finally, plot the frequency response of the output samples. As expected, we essentially only get the peak at 730 Hz
Y = fft(yval)*Ts; fy = (0:numel(Y)-1)/numel(Y)/Ts;
figure
plot(fy,abs(Y)),xlim([0 1000]),xlabel('Hz');ylabel('Amplitude')
The preceding analysis focused only on amplitude. Further adjustments would need to be made in the analysis if we want Y to have correct phase for comparison to the phase of the CTFT of y(t), i.e., Y(w).
  6 commentaires
Image Analyst
Image Analyst le 28 Juil 2023
Modifié(e) : Image Analyst le 28 Juil 2023
@Paul the issue is he's not using MATLAB. He's using Octave - an open source MATLAB look-alike. At least he's tagged the post as Octave. But Octave is not the same as MATLAB in all respects. He should ask in an Octave forum, not here.
Paul
Paul le 29 Juil 2023
I don't think that's the issue at all, at least as far as why he didn't get the result he was expecting.
Also, the question did explicitly say "and plotted it in matlab."
If a question has an octave tag, which I didn't see because by and large I don't look at the tags, should it not be answered and perhaps closed instead?

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Octave dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by