# Can you help me how to move forward

2 views (last 30 days)
kraken on 21 May 2021
Edited: William Rose on 21 May 2021
Fs = 5000;
Fn = Fs/2;
numSamples = 250*Fs;
t = linspace(0,50,numSamples);
x = sin(2*pi*100*t) + 0.7*cos(2*pi*500*t);
xn = x + 5.*randn(size(t));
originalSpectrum = fft(x);
noisySpectrum = fft(xn);
Fv = linspace(0, 1, fix(numel(t)/2)+1)*Fn;
Iv = 1:numel(Fv);
figure
subplot(4, 1, 1);
plot(Fv, abs(originalSpectrum(Iv))*2, 'b-', 'LineWidth', 2);
grid on;
xlabel('t');
ylabel('x');
subplot(4, 1, 2);
plot(Fv, abs(noisySpectrum(Iv))*2, 'b-', 'LineWidth', 2);
dee = designfilt('lowpassfir', 'PassbandFrequency', 500, ...
'StopbandFrequency', 600, ...
'PassbandRipple', 1, 'StopbandAttenuation', 100, ...
'SampleRate', 5000);
dataOut = filter(dee,xn);
subplot(4, 1, 3);
plot(dataOut);
trfcn = fft(dataOut)./noisySpectrum;
subplot(4,1,4)
plot(Fv, abs(trfcn(Iv)))
title('Transfer Function Of Filter')
frequencies are transferred to the output at different speeds, there is a time difference
between input and output signals. A method to compensate for this time difference
It is wanted to be followed as such.
1-Calculate the mean value of the group delay of the designed LPF in Matlab. Let us define this value as 𝐷.
2- Before applying the noisy signal to the LPF input, do zero-padding on this signal for a number of samples equal to 𝐷.
3- Apply this signal to the LPF input.
4-Discard the first D samples of the signal at the LPF output to obtain the denoised and delay compensated signal 𝑦𝐿𝑃𝐹−𝐶 [𝑛].
Can you help me how to move forward

William Rose on 21 May 2021
Edited: William Rose on 21 May 2021
@krakenThe code you provided runs - good.
Are you aware that, because of how you constructed t, the sampling rate in t does not equal Fs? The sampling rte of t is 25 kHz. You can confirm that by inspecting the values in t. Is that your intention? Therfore, when you construct the digital filter, and you pass sampling rate of 5000, that is not the sampling rate of x. That is the first thing I would correct.
The reason this happens is that you define
numSamples = 250*Fs;
which a number corresponding to 250 seconds worth of samples at rate Fs (which is 5000). Then you do
t = linspace(0,50,numSamples);
which creates a vector t that goes from 0 to 50 seconds, and has 250*5000 elements. Therefore the sampling rate of t is 50/1249999=24.99998 kHz. I don't use linspace() to make time vectors, because if you want a normal sampling rate such as 5 kHz, and if you want an even number of samples (which is nice for FFTs), and if your first time value is zero, then you do not want the last sample to be exactly a round number such as 50 seconds. You want the last sample to be one less than the round number. So, if you are sampling at 5 kHz, then you want the last sample to be 49.9998 seconds.
I don't know what you intended: 250 seconds of data at 5 kHz, or 50 seconds at 25 kHz, or 50 seconds at 5 kHz. I will assume the last one is what you wanted. Therefore I recommend changing lines 3,4 of your code to the following.
numSamples = 50*Fs;
t = (0:numSamples-1)/Fs;
That's a start.

William Rose on 21 May 2021
@kraken, next part of my answer:
ylabel('x');
xlabel('t');
but the correct labels should be
ylabel('|X(f)|');
xlabel('Frequency (Hz)');
Your plot 3 is the time domain noisy signal after filtering, but you do not use the time vector when plotting. I recommend that you do, so the viewer knows the true time scale. Also, include axis labels for the same reason.
%plot 3: noisy signal after filtering
subplot(4, 1, 3);
plot(t,dataOut);
xlabel('Time (s)');
ylabel('y(t)');
With the changes so far, I now have the following graphical output:

William Rose on 21 May 2021
gdel=grpdelay(dee);
returns the group delay at 512 points from frequency 0 to π rad/sec. dee has order 159, which you know from doing
>> disp(filtord(dee))
159
The Matlab help recommends that one evaluate the delay at a number of freuqncies greater than the filter order, and 512>159 so you're OK there. It turns out the delay is the same at all frequncies - which is typical for symmetric FIR filters, and the mean, rounded to the nearest integer, is
>> D=round(mean(grpdelay(dee)))
D =
80
This corresponds to a delay of D/Fs=80/5000=0.016 sec = 16 msec. We will zero-pad xn(t) at the end, not the beginning, because if we zero-pad the beginning, there will be even more delay. So we add zeros at the end of xn, then filter, then discard the first D samples from the fltered signal.
dataOut = filter(dee,[xn,zeros(1,D)]);
dataOutDC = dataOut(D+1:end); %delay-compensated version of the filter output
Then plot. I made three plots of the time domain signals xn(t), dataOut(t), and dataOutDC(t), with increasing expansion of the x-axis limits. IN the last plot you can see the initial delay of the filtered output, before delay complensaiton, and you can see that the delay-ompensated signal fits nicely on top of the unfiltered noisy signal.
By the way, the way you have computed the transfer function estimate, , for the filter, from the input and output data, is a good first guess for how to do it, but it is not the recommended way. Your estimate of the TF is rather noisy. A much preferable way is
where * indicates complex conjugate. There are books devoted to this subject of transfer function esitmation, which explain why the formula above is better.
William Rose on 21 May 2021
Here is a plot of the filter transfer function, estimated from the input and output, using the method which I mentioned above:
which you can implement using cpsd() in the SigProc toolbox (here):
>> gest=cpsd(xn,dataOut(1:end-D))./cpsd(xn,xn);
>> f=linspace(0,Fn,length(gest));
>> plot(f,abs(gest),'r-');
This estimate of the transfer function is much less noisy than your estimate, and it appears consistent with the filter design specifications.