1
$\begingroup$

I'm trying to plot some sinusoidal signals in Matlab.

But while frequency is getting higher (closing to fs/2), results are getting distorted.

I guess it's lack of my knowledge but distortion is starting soon. My peak points are getting away from the original amplitude.

What is the reason and how can I correct it?

Code is;

fs =441; %discrete *100 - Hz
%data
phase = 0;
A = 6;
f = 50; % *100 -Hz
%sampling
n=0:2/f*fs;
Ts=1/fs;
Tn=n*Ts;
%signal
Xn=A*cos(2*pi*f*Tn+phase*pi);
%plotting
stem(Tn,Xn)
figure;plot(Tn,Xn)

Normal signal; 1kHz @44100Hz enter image description here

Distorted signal #1; 5kHz @44100Hz enter image description here

Distorted signal #2; 15kHz @44100Hz enter image description here

Distorted signal #3; 20kHz @44100Hz enter image description here

Almost ok signal; 22050Hz @44100Hz enter image description here

$\endgroup$
1

2 Answers 2

2
$\begingroup$

One way to avoid low time resolution is to upsample signals before whatever processing.

To answer this question I have modified the Mathworks example that downsamples DAB (audio) signals to match DVD (audio) rates, example available here:

https://uk.mathworks.com/help/signal/ref/upfirdn.html

So Let's get started

close all;clear all;clc

1.- Signal parameters

fs1=44.2e3;Ts1=1/fs1;     % 22.624us 
f1=5e3;T1=1/f1;         % .2ms  

spf=500;  % samples per frame

2.- Low pass filter parameters

f_pass=21e3;f_stop=22e3;  % LPF start and stop frequencies have to be < fs1
Rp=1; % [dB] pass band ripple
Lstop=60; % [dB] pass band attenuation

3.- Upsampling parameters

fs2=4*fs1  % upsampling frequency
[L,M] = rat(fs2/fs1)
f_ = (fs1/2)*min(1/L,1/M)

3.- Generating signal s1 same as Sine1()

3.1.- s1 values

Sine1 = dsp.SineWave('Frequency',f1,'SampleRate',fs1,'SamplesPerFrame',spf);
s1=Sine1();

3.2.- s1 time reference

t1=[0:1:spf-1]*Ts1;

4.- Building low pass filers

4.1.- LPF1 with command dsp.LowpassFilter

LPF1= dsp.LowpassFilter('PassbandFrequency',f_pass, ...
                                     'StopbandFrequency',f_stop, ...
                                     'PassbandRipple',Rp , ...
                                     'StopbandAttenuation',Lstop, ...
                                     'SampleRate',fs2,  ...
                                     'FilterType','IIR');

fvtool(LPF1,'Analysis','freq')  % check H(f) of LPF1

enter image description here

4.2.- LPF11 another low pass filter with quite the same parameters as LPF1 and built with alternative command designfilt

% generating anti-aliasing filter LPF11
 LPF11= designfilt('lowpassfir', ...
                         'PassbandFrequency',.9*f_, ...
                         'StopbandFrequency',1.1*f_, ...
                         'PassbandRipple',3, ...     % Rp, ...
                         'StopbandAttenuation',40, ... % Lstop, ...
                         'DesignMethod','kaiserwin', ...
                         'SampleRate',fs2);
fvtool(LPF11,'Analysis','freq')  % check H(f) of LPF11

enter image description here

5.- Building Spectrum Analyzers

SA = dsp.SpectrumAnalyzer('PlotAsTwoSidedSpectrum',false, ...
                                          'SampleRate',4*fs1, ... % Sine1.SampleRate, ...
                                          'FrequencyResolutionMethod','WindowLength', ... 
                                          'FrequencyVectorSource','Property',...
                                          'FrequencyVector',[-5000 000],...
                                          'NumInputPorts',2,...
                                          'FrequencyResolutionMethod','RBW', ...
                                          'RBW',40, ...
                                          'FFTLength',4096, ...
                                          'ShowLegend',true, ...
                                          'YLimits',[-30,45]);

% an alternative Spectrum Analyzer I used to do some checks
% SA2 = dsp.SpectrumAnalyzer('PlotAsTwoSidedSpectrum',false, ...
%                                           'SampleRate',4*fs1, ... % Sine1.SampleRate, ...
%                                           'NumInputPorts',2,...
%                                           'ShowLegend',true, ...
%                                           'FrequencyResolutionMethod','WindowLength', ...  % 'RBW'
%                                           'FFTLengthSource','Auto', ... % 'Property', ...
%                                           'RBW',20, ...
%                                           'FFTLength',2048, ...
%                                           'WindowLength',128, ...
%                                           'TimeResolution',1e-4, ...
%                                           'YLimits',[-30,45]);

6.- Filtering s1

6.1.- Without upsampling : s1 through filter LPF1

I have added a small amount of nose to s1 and then it's useful to do some averaging

SA.ChannelNames = {'s1','LPF(s1)'};
y1avg=zeros(1,numel(s1));
for i = 1 : spf
s1n=s1; % + 0.1.*randn(Sine1.SamplesPerFrame,1);
    y1 = LPF1(s1n);
    SA(s1n,y1); % filtering 
     y1avg=.5*(y1avg+y1); % averaging  
end
release(SA)

enter image description here

figure
plot(t1,y1avg,'b');
hold on;grid on
stem(t1,y1avg,'b');
plot(t1,s1','r');
stem(t1,s1,'r');
xlim([10*T1 13*T1]);
xlabel('t1');title('s1(t) y2=avg(LPF(s1)) no ups-LPF1')

enter image description here

6.2.- With upsampling : s1 through filter LPF1

h0=L*tf(LPF11);
s1up = upfirdn(s1,h0,L,M);
d1 = floor(((filtord(LPF11)-1)/2-(L-1))/L);
s1up = s1up(d1+1:end);
t1up = (0:(length(s1up)-1))/fs2;

figure
stem(t1up,s1up,'*')
hold on;grid on;
stem(t1,s1);
xlim([10*T1 13*T1]);
legend('s1 original','s1 resampled','Location','southeast')
hold off

enter image description here

upsampled s1 before filter

SA.ChannelNames = {'s1','LPF(s1)'};
y1upavg=zeros(numel(s1up),1);
for i = 1 : spf
    y1up = LPF1(s1up);
    y1upavg=[y1upavg y1up];
    SA(s1up,y1up);
end
y1upavg(:,1)=[];  
y2up=mean(y1upavg,2); % averaging
release(SA)

enter image description here

original s1 and averaged upsampled antialiased y2 .

Antialiasing is another way to say low pass filtered when the low pass filter is trimmed to the input signal.

figure
plot(t1up,y2up,'b');
hold on;grid on
stem(t1up,y2up,'b');
plot(t1,s1','r');
stem(t1,s1,'r');
xlim([10*T1 13*T1]);
xlabel('t1');title('s1(t) y2=avg(LPF(s1))')

enter image description here

Repeating for s2 s3 s4 shows that s1 frequency can be really close to a good 'flat' filter frequency edges and yet avoid amplitude degradation if the signal is upsampled just before 1st filter.

In analog signal processing, as example RF microwave filtering with transmission-line built filters there's no need to 'condition' signals as shown here with upsamping, but it's certainly important to use anti-aliasing filters nevertheless.

I have included LPF11 in addition to LPF1 because while LPF1 has the question specs and a cut-off half of the sampling frequency fs1, LPF11 has cut-off just above the input carrier s1 ; which is the best way to perform low pass filtering.

There's a lot of spectrum between f1 and fs1/2 22kHz that is good practice to avoid into a low pass filter intended for carrier f1 at just 5kHz .

Thanks for reading

$\endgroup$
4
$\begingroup$

It's simply a drawing artifact. You are plotting the wave form by connecting the individual samples with straight lines. That's wrong.

The values between samples can be determined through sinc interpolation, not through linear interpolation. The difference between the two becomes more and more visible as the frequency increases.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.