Can you help me how to move forward

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

Answer 1

[William Rose]

@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.

Answer 2

[William Rose]

@kraken, next part of my answer:

Your plot 1 uses

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:

Answer 3

[William Rose]

@kraken,

Now for the answers to your specific quesitons:

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]

@kraken, your code includes

plot(Fv, abs(originalSpectrum(Iv))*2, 'b-', 'LineWidth', 2);

Why do you have a factor of 2? Because you are plotting the one-sided spectrum? If so, then multiplying all the amplitudes by 2 is not the way to go, for a couple of reasons.

• You don't multiply the first (f=0) and last (f=fNyquist) by two, because they are not "doubled". (If the numSamples is odd, the last element of the spectrum will not be at the Nyquist frequency, so it should be doubled in the one-sided power spectrum. However, numSamples is even in your case, and in the great majority of cases, so the last element of the spectrum is at the Nyquist freq, and therefore should not be doubled in the power spectrum.)
• The power, not the amplitude, should be multiplied by two.

[William Rose]

@kraken,

Are you a fan of giant sea monsters? Or Seattle hockey?

Code to make the figures in my answer is attached.

[William Rose]

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.