Verfasst am: 11.11.2014, 14:31
Titel: Upsampling mit anschließendem Tiefpass
Hallo zusammen,
derzeit versuche ich, die Interpolationsfunktion von Matlab (interp(..)) für meinen Zweck selbst nachzubauen, um mir ein besseres Verständnis zu erarbeiten.
Hierzu habe ich ein Skript File entwickelt, in dem ich Schritt für Schritt das Upsampling durchführe.
Einmal lasse ich das Skript File wie folgt ablaufen:
Code:
x = [167910142795392];
upsamplevalue = 1;
Fs = upsamplevalue; % Sampling frequency
x = upsample(x, upsamplevalue);
%x = interp(x,2);
T = 1/Fs; % Sample time
L = length(x); % Length of signal
t = (0:L-1)*T; % Time vector
% f_tp = 0.5; % %Tiefpassfilterung % [b a] = butter(50,f_tp/(0.5*Fs),'low'); % Koeffizienten der Übertragungsfunktion, f_tp ist die Grenzfrequenz, fs die Abtastfrequ., Butterworthfilter % figure(3) % freqz(b,a); % x = 2*filter(b,a,x); %x ist das Signal
% h = fir1(100,0.5); % % x = conv(x,h);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
X = fft(x,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
Dabei erhalte ich zwei Plots. Einmal meine Abtastwerte im Zeitbereich, und einmal das dazugehörige Spektrum (1.Nyquistzone).
Das Spektrum startet dabei mit einer maximalen Amplitude von ca. 25.
Anschließend lasse ich das Skript erneut durchlaufen, erhöhe aber meinen Parameter upsamplevalue auf den Wert 2.
Dadurch habe ich zwischen allen Abtastwerten zusätzliche Nullstellen. Das Spektrum ist dabei sozusagen zweimal innerhalb meiner 1.Nyquistzone zu sehen, da diese aufgrund des Upsamplings und der folgich doppelt so großen Abtastrate verdoppelt wurde.
Hier mein erstes Problem, mein Spektrum sollte meinen Theorieerkenntnissen nun nur doppelt vorhanden sein, jedoch mit gleicher Amplitude. Bei mir ist die Amplitude aber nun nur noch betragsmäßig halb so groß (maximal Wert zu Beginn bei ca. 13)
Die FFT Funktion habe ich von Matlab beispielen verwendet, dabei wird die fft wie im Code zu sehen ist, durch L geteilt (Länge des Signalvektors).
Warum ist mein Spektrum beträgsmäßig halb so groß nach dem Zerostuffing, Upsampling?
Zudem wollte ich nun als letzten Schritt, mein upgesampeltes Signal noch mit einem Tiefpass filtern. Hierzu kommentiere ich aus meinem Quellcode zusätzlich die beiden Zeilen ein, zur Erzeugung des fir1 Tiefpasses mit der Grenzfrequenz 0.5 und falte den Tiefpass anschließend mit meinem Eingangssignal x.
Im Spektrum erkenne ich dabei auch, dass die Frequenzen unterhalb von 0.5 Hz weggeschnitten wurden.
Jedoch wurden in meinem Zeitsignal wie erwartet nicht nur die Nullstellen interpoliert sondern auch meine ursprünglichen Stützstellen verändert.
Woran liegt das, verwende ich das Filter falsch?
Wenn ich die von Matlab verwendete interp-Funktion verwende erhalte ich auch das erwartete Spektrum, aber die korrekte Zeitfunktion.
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
X = fft(x,NFFT);
f = Fs/2*linspace(0,1,NFFT/2+1);
% Berechnung des Amplitudengangs aus dem komplexen Frequenzvektor H:
amplH = abs(X);
% Darstellung des interessierenden Frequenzbereichs des % Amplitudengangs (0...fn) und % daran angepasste Amplitudenskalierung (Normierung auf N/2):
amplitudengang = [amplH(1)/NFFT amplH(2:NFFT/2)/(NFFT/2) amplH(NFFT/2 + 1)/NFFT]; % DC-Bin auf N normieren!
Ich habe auch unter dem angegebenen Link nochmals nachgeschaut, und auch deine Änderungen eingepflegt, aber weiterhin ist das Problem, dass bei einem upsamplevalue von 2 die gesamte Amplitude um den Faktor 2 kleiner ist. (Also im Vergleich zu einem upsamplevalue von 1, also kein upsampeln).
Ich könnten natürlich auch einfach immer auf meinen Maximalwert der FFT normieren, jedoch würde ich dann nur die Amplitude "schmutzig" gleichsetzen, obwohl dies gar nicht stimmt.
Warum ist also meine Amplitude um den Faktor 2 skaliert, bei einem upsamplevalue von 2? In der Theorie sollte doch das gleiche Spektrum erscheinen, nur dass die Nyquistzone doppelt so groß ist, und ich somit ein Spiegelspektrum erhalte.
Hat desweiteren noch jemand eine Idee, wie ich das zweite Spiegelspektrum z.B. mit einem FIR Filter wegbekomme. Ich habe dies bislang versucht, jedoch komme ich auf keinen grünen Zweig.
Nicht sehr schön...aber dafür die Lösung. Ich habe mal als Testsignal einen einfachen Sinus erzeugt, da mir damit der Zusammenhang klarer ist. Das Signal muss um den Faktor upsamplevalue vor dem Filter multipliziert werden, da ja eine lineare Interpolation gefordert wird...sprich der Mittelwert zwischen zwei real existierenden Messwerten gesucht wird. Durch das Einfügen von Nullen, wird ja aber der Mittelwert kleiner und der Korrekturfaktor somit unumgänglich.
upsamplevalue = 2;
Fs = 100; % Sampling frequency
L = 1024;
T = 1/Fs; % Sample time
t = (0:L-1)*T; % Time vector
df = Fs/L; % Auflösung Freq.-spektrum
freq = df * 100; % ganzes Vielfaches von df damit kein leakage auftritt
x = 10 * sin(2*pi*freq*t); % sinus mit Freq. = 10.24Hz und Amplitude = 10
Fs_up = Fs * upsamplevalue;
x_up = upsample(x, upsamplevalue); % enspricht Einfügen von upsamplevalue Nullen zwischen 2 Werten
%x = interp(x,2);
T_up = 1/Fs_up;
L_up = length(x_up); % Length of signal
t_up = (0:L_up-1)*T_up; % Time vector
df_up = Fs_up/L_up; % Auflösung Freq.-spektrum
y_interp = interp(x,upsamplevalue);
f_tp = 0.5*Fs;
% %Tiefpassfilterung [b, a] = butter(4,f_tp/(0.5*Fs_up),'low'); % Koeffizienten der Übertragungsfunktion, f_tp ist die Grenzfrequenz, fs die Abtastfrequ., Butterworthfilter
%b = fir1(10,f_tp/(0.5*Fs_up));
% figure(3) % freqz(b,a);
x_up_filt = filtfilt(b,a,upsamplevalue*x_up); %x_up ist das Signal nach upsample
% filtfilt wird verwendet damit kein Filterdelay zu sehen ist % Korrekturfaktor = upsamplevalue
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
X = fft(x,NFFT);
f = Fs/2*linspace(0,1,NFFT/2+1);
% Berechnung des Amplitudengangs aus dem komplexen Frequenzvektor H:
amplH = abs(X);
% Darstellung des interessierenden Frequenzbereichs des % Amplitudengangs (0...fn) und % daran angepasste Amplitudenskalierung (Normierung auf N/2):
amplitudengang = [amplH(1)/NFFT amplH(2:NFFT/2)/(NFFT/2) amplH(NFFT/2 + 1)/NFFT]; % DC-Bin auf N normieren!
subplot(2,1,2);
NFFT = 2^nextpow2(L_up); % Next power of 2 from length of y
X = fft(x_up_filt,NFFT);
f = Fs_up/2*linspace(0,1,NFFT/2+1);
% Berechnung des Amplitudengangs aus dem komplexen Frequenzvektor H:
amplH = abs(X);
% Darstellung des interessierenden Frequenzbereichs des % Amplitudengangs (0...fn) und % daran angepasste Amplitudenskalierung (Normierung auf N/2):
amplitudengang = [amplH(1)/NFFT amplH(2:NFFT/2)/(NFFT/2) amplH(NFFT/2 + 1)/NFFT]; % DC-Bin auf N normieren! stem(f,amplitudengang,'r'); hold on;
X = fft(y_interp,NFFT);
% Berechnung des Amplitudengangs aus dem komplexen Frequenzvektor H:
amplH = abs(X);
% Darstellung des interessierenden Frequenzbereichs des % Amplitudengangs (0...fn) und % daran angepasste Amplitudenskalierung (Normierung auf N/2):
amplitudengang = [amplH(1)/NFFT amplH(2:NFFT/2)/(NFFT/2) amplH(NFFT/2 + 1)/NFFT]; % DC-Bin auf N normieren!
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS