# Nyquist frequency Plotting Distortions

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

Distorted signal #1; 5kHz @44100Hz

Distorted signal #2; 15kHz @44100Hz

Distorted signal #3; 20kHz @44100Hz

Almost ok signal; 22050Hz @44100Hz

• Related Oct 18, 2022 at 10:37

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


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


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)


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')


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


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)


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))')


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 .