2 views (last 30 days)

Show older comments

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:

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:

William Rose
on 21 May 2021

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
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:

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

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

Start Hunting!