Verfasst am: 29.05.2012, 17:49
Titel: Cepstrum wieder in Spektrum umwandeln
Hallo,
ich sitze gerade an einer Funktion, die ein Audiosignal einliest, das Cepstrum berechnet und an dem dann einige Veränderungen vorgenommen werden sollen. Als Ausgabe soll wieder ein Audiosignal erstellt werden.
Zum Testen habe ich das Cepstrum erstmal unverändert übergeben.
Zur Überprüfung lasse ich den Amplitudengang des eingelesenen und des aus dem Cepstrum rückgerechneten Signals plotten. Die Plots liefern augenscheinlich dasselbe Ergebnis.
Aber ein Audiosignal kann ich daraus nicht wiederherstellen. Ersetze ich im Code zur Erstellung der Wavedatei 'neu_fft' durch 'x_fft', funktioniert das jedoch.
Daher vermute ich, dass bei der Umwandlung des Cepstrums zurück in das Spektrum etwas nicht stimmt.
Für eure Anregungen und Hilfen wäre ich dankbar!
Code:
function[] = test() % Signal einlesen [signal fs] = wavread('./audio/test.wav');
% Parameter
N = 1024;
freq = 0:fs/N:fs/2;
% Spektrum und log-Amplitudengang des Signals
x_fft = fft(signal);
x_amp = abs(x_fft);
x_amp_scale = [x_amp(1)/N; x_amp(2:N/2)/(N/2); x_amp((N/2)+1)/N];
x_log = 20*log10(x_amp_scale);
% Cepstrum berechnen
signal_ceps = rceps(signal);
% sollte gleichbedeutend sein mit: % signal_ceps = real(ifft(log(abs(fft(signal)))));
nachdem die Wiederherstellung des Audiosignals nun ja funktioniert, möchte ich gerne mit dem Cepstrum arbeiten.
Ich möchte das Cepstrum aus zwei verschiedenen zusammensetzen:
Den einen Teil aus Signal X, den anderen aus Signal Y. Ein threshold bestimmt die Teilung.
Wenn ich DANN allerdings das Audiosignal wiederherstellen will, wird das ganz seltsam:
nehme ich nur den ersten Koeffizienten von Signal X, ist der output derselbe wie Y, jedoch wesentlich leiser
für alle anderen thresholds ist der Output nicht abspielbar, und das Spektrum liegt in vierstelligen Amplitudenbereichen (in dB!).
Hat jemand eine Idee, woran das liegen könnte? Falsche Skalierung?
Mit icceps bekomme ich keine vernünftigen Ergebnisse.
Zur Erstellung des Cepstrums nutze ich mittlerweile nicht mehr rceps(), sondern ifft(log(abs(fft(signal)))).
Die Wiederherstellung erfolgt dann per exp(fft(out_ceps,fft_size)) und liefert offensichtlich das richtige Spektrum.
Ich muss zugeben, dass dieses Thema für mich Neuland ist. Ich lese mich da gerade erst ein. Gleichzeitig musste ich aber auch feststellen, dass diese Herangehensweise mir nicht wirklich weiterhilft. Ich wollte mit dem Cepstrum mein Spektrum glätten. Da meine Signale aber keine signifikanten Stellen im Cepstrum aufweisen (es handelt sich nicht um ein Audiosignal), nutzt mir der Vorteil, welcher durch die Logarithmierung entsteht leider nichts. Oder ich übersehe hier auch noch einiges...
Nun aber zu meinem Einwand. Ich muss dir Recht geben, dass die Funktion icceps leider nicht das richtige Ergebnis liefert, wenn man dieses zurückgewandelte Spektrum wieder in den Zeitbereich transformiert. So weit ich das nun verstanden habe, kommt aber ohnehin nicht das identische Zeitsignal wieder, da ja die Phaseninformation durch das Transformieren ins Cepstrum verloren geht. Allerdings bin ich doch sehr verwundert über des Betragsspektrum von s2_lift2 nach der Rücktransformation aus dem Cepstrum. s2_lift hingegen, stimmt mit dem original Spektrum überein. Wandelt man dieses s2_lift aber wieder in den Zeitbereich entsteht "Müll", während s2_lift2 bis auf eine Phasenverschiebung mit s2 übereinstimmt.
Hier mein code...als Testsignal habe ich einfach ein sinus mit Echo = s2. Die Funktion FFT_h_da hast du zwar nicht, da passiert aber nichts weiter aufregendes. Es ist mit dem "umfassenden FFT Bsp" hier in der Skriptecke zu vergleichen.
Code:
clear;
Ts = 0.001;
Fs = 1/Ts;
N = 1024;
df = Fs/N;
t = 0:Ts:(N-1)*Ts;
f = 20*df;
s1 = sin(2*pi*f*t);
% Echo
s2 = s1 + 0.5*[zeros(1,256) s1(1:768)];
% Lifter % Test mit Hamming Fenster
indx = find(t<0.235,1,'last');
win = rectwin(2*indx);
%win = hamming(2*indx);
lifter = [win(indx:end); zeros((N/2)-indx-1,1)];
clifter = cceps(lifter);
c_lift = c;%(c(1:N/2) .* lifter');
% Glätten des Spektrums durch "Tiefpassfilterung des Cepstrums" % einfacher herangehensweise durch Löschen der hohen Anteile = % Rechteckfenster
%c_lift(indx+1:end-indx+1) = 0;
Auf den beiden Plots sieht man aber auf was ich hinaus will. Rot und blau stimmen im Freq.bereich überein. Das grüne Spektrum kann ich aber überhaupt nicht nachvollziehen. Bei der Rücktransf. in den Zeitbereich liefert dann aber dieses merkwürdige Spektrum annähernd s2 wieder, während aus dem roten Spektrum "Müll" entsteht.
Also ich habe gestern auch nochmal hin und her probiert mit den verschiedensten Möglichkeiten der Implementierung für das Cepstrum. Das Zusammenstellen und Rücktransformieren funktioniert jedoch mittlerweile. Allerdings nicht sehr zufriedenstellend, weil das Rauschen immernoch präsent ist, da fällt mir nur gerade nichts zu ein.
Glätten und Audiosignale sind ansonsten auch mein Arbeitsbereich
Im Moment glätte ich mit dem in Matlab integrierten LPC, davor hatte ich mich aber auch anch Cepstrum-Glättung umgesehen und dazu eine ganz hilfreiche Masterarbeit gefunden:
www.rs-met.com/documents/dsp/Magisterarbeit.pdf
Den darin beschriebenen Algorithmus in der verbesserten Form habe ich mal nachgebaut und der lief auch sehr ordentlich. Ich hab ihn dir mal angehangen.
Wenn du mir den Code der fft_h_da Funktion zeigst, kann ich deinen Code auch laufen lassen.
Danke für den Link und den m-file...ich werde mir das später mal genauer anschauen.
Die Funktion FFT_h_da kann ich nicht veröffentlichen, da hier noch einige andere Sachen über das Betragsspektrum hinaus erfolgen. Ich habe mal nur zur Darstellung des Spektrums die Fkt. verkürzt.
Code:
function[mag, mag_dB, fv] = FFT_betragsspektrum( signal, nfft, fa) ;
% Input: % signal = Signal im Zeitbereich % nfft = Anzahl Messwerte für fft % wenn nfft > length(sig) -> fft(sig,nfft) führt Zeropadding durch % fa = Abtastfreq. % Output: % Magnitude des Spektrums linear und dB skaliert % Frequenzvektor fv in [Hz] von 0...fa/2
% un-,gerade Anzahl Messwerte? ifmod(nfft,2) == 0;
k = (nfft/2) + 1;
else
nfft = nfft + 1;
k = (nfft/2) + 1;
end
Die Masterarbeit hat mir ehrlich gesagt überhaupt nicht weitergeholfen.
Hast du hier noch neue Erkenntnisse sammeln können bzw. kannst mir evtl. meine Fragen beantworten? Das rücktransformierte Spektrum aus dem Cepstrum ergibt für mich keinen Sinn (siehe Plots). Ebenso die Ergebnisse der Rücktransformation in den Zeitbereich.
Kann man denn überhaupt das Zeitsignal nach den Transformationen wieder herstellen?
pullmoll89
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 16.07.2012, 11:10
Titel:
Hey.
Ja, ich habe mittlerweile herausgefunden, wie es funktioniert. Zumindest läuft das bei mir so richtig^^
Um die Funktion cceps() zu beutzen musst du folgende folgende Änderungen machen:
Die Spektren von s2 und s2_lift sehen dann immer noch so aus wie in deinem Plot, d.h. sie sind gleich.
s2_lift ist dann aber schon das Zeitsignal, muss also nicht weiter behandelt werden [ real(ifft(...)) ] und stimmt dann im letzten Plot mit s2 überein.
Ich berechne meine Glättung bzw. Einhüllende allerdings mit der exp()-Variante.
Das Cepstrum meines Ausgangssignals füge ich dann aus Einzelteilen zusammen, wobei du beachten musst, dass das Cepstrum zweiseitig ist, d.h. die niedrigen Koeffizienten befinden sich z.b. links und rechts (wie in der Berechnung oben).
Und mit
Ich habe gestern abend mal deine Vorschläge getestet.
Mit der Erweiterung für cceps bzw. icceps stimmen nun die Ergebnisse.
Allerdings erreiche ich mit den Einzelschritten nicht wieder das orig. Signal im Zeitbereich nach den Transformationen.
Ich nehme mal an, dass es mit der Phaseninformation zu tun hat, die du hier ja gar nicht berücksichtigst. Offensichtlich bist du daran ja auch gar nicht interessiert. In der Doku zu cceps wird ja auch der interne Weg der Funktion beschrieben, wobei eine Funktion rcunwrap verwendet wird, die in MatLab gar nicht zugänglich ist.
In so fern wundert es nicht, dass ich das Zeitsignal wie in meinem Codebsp. nicht wieder rekonstruieren kann.
Code:
clear;
Ts = 0.001;
Fs = 1/Ts;
N = 1024;
df = Fs/N;
t = 0:Ts:(N-1)*Ts;
f = 20*df;
s1 = sin(2*pi*f*t);
% Echo
s2 = s1 + 0.5*[zeros(1,256) s1(1:768)];
figure(1) plot(t,s1,t,s2) xlabel('Zeit in s') ylabel('Amplitude') legend('original Signal ohne Echo','original Signal mit Echo') title('Zeitbereich') grid on;
figure(2) % Darstellung des Cepstrums plot(1:length(signal_ceps),20*log10(signal_ceps),'b',1:length(signal_ceps_low),20*log10(signal_ceps_low),'r--');
xlabel('Quefrency');
ylabel('amplitude cepstrum') title('cepstrum') legend('Cepstrum von s2', 'signal_ceps_low') grid on;
figure(3) % Spektrum von s2 darstellen [mag, mag_dB, fv] = FFT_betragsspektrum( s2, N, Fs) ;
plot(fv,mag_dB,'b.-') ; hold on;
% Spektrum der Signale zurück aus Cepstrum plot(fv,signal_env,'r--') ; hold on;
plot(fv,env_log,'g--') ; hold on;
%axis([0 Fs/2-12040]);
xlabel('Frequenz in [Hz]') ;
ylabel('Magnitude in [dB]') ;
legend('original Signal mit Echo','Signal mit Rechteckfenster bearbeitet','Env. Funktion') title('Betragsspektrum') grid on;
figure(4) % Wieder in Zeitbereich - im Cepstrum verändertes Signal
s2_lift = ifft(signal_fft_low,fft_size);
plot(t,s1,'g',t,s2,'b',t(1:length(s2_lift)),s2_lift,'r.-');
xlabel('Zeit') ylabel('Amplitude') legend('original Signal ohne Echo','original Signal mit Echo','Signal - Cepstrum bearbeitet') title('Zeitbereich') grid on;
Ich habe auch mal die Funktion Envelope getestet. Wobei ich damit auch nichts anfangen kann. Wie ich schon schrieb, befasse ich mich nicht mit Audiosignalen und benötige deshalb auch keine Einhüllende meines Spektrums. Ich bin an einer Glättung meines Spektrums interessiert. Durch das einfache Löschen von Bereichen, wie du es vornimmst, komme ich aber auch nicht weiter. Außerdem fehlt mir ein Kriterium für die Wahl des Threshold. Hier wird meine Signalfreq. von s1 eher noch verstärkt, was ich nicht gebrauchen kann. Mir geht es um die Trennung/ Verarbeitung von Signalen, die multiplikativ verbunden sind. Durch den Wechsel ins Cepstrum wird dieser Zusammenhang ja zu einer Addition und somit ist eine Filterung wiederum möglich. Ich benötige hier aber andere Filtertypen und konnte leider für das Lifting wie auch die Erstellung des Lifters keine Bsp. finden, weshalb ich nun danach gefragt habe.
Erstellt man hier einfach die Impulsantwort des Filters im Zeitbereich und transformiert sie dann ins Cepstrum?
Ich hoffe du kannst mir hier noch weiterhelfen...aber dennoch vielen Dank für die bisherige Unterstützung.
pullmoll89
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 17.07.2012, 15:05
Titel:
Also was den Lifter angeht, da kenne ich mich gar nicht mit aus.
Aber ich glaube, du kannst das Signal generell nicht wiederherstellen, wenn du nur die signal_ceps_low nimmst, und den Rest zu Null setzt.
In meinem Falle berechne ich so die Einhüllende ohne Probleme, aber zur Rekonstruktion tausche ich die unteren Cepstrum-Koeffizienten gegen andere aus, sodass alle wieder da sind.
Dann habe ich bei mir auch keine Probleme mit der Rekosntruktion im Zeitbereich.
In diesem Teil löscht du in der Mitte des Signals ab der Stelle threshold. Da der Vektor wie bei der FFT für ein beidseitiges Cepstrum (Spektrum) gilt, ist der Vektor gespiegelt und das Ende entspricht dem Anfgang. Deshalb bleibt das Ende ab end-threshold auch unverändert.
Da hier aber nur der Betrag von signal_fft für signal_ceps genutzt wird, kann die Phaseninformation nicht wieder hergestellt werden. Was du hier nicht zeigst, wie out_ceps zusammengestellt ist.
rekonstruiert das Zeitsignal nicht mehr, selbst wenn Threshold so gewählt wird, dass signal_ceps_low = signal_ceps ist und damit das Signal unverändert bleibt.
pullmoll89
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 17.07.2012, 20:07
Titel:
Mit dem signal_ceps_low berechne ich nur die Einhüllende. Dazu genügt das zu-Null-setzen der höheren Koeffizienten, die sich (da zweisietig) ja in der Mitte befinden. Höhere Koeffizienten im Cepstrum bedeuten ein feineres Signal mit mehr Artefakten etc.
"Austausch" heißt in meinem Falle, dass ich die niedrigen Koeffizienten ersetze. Und zwar habe ich als Input ein Sprachsignal mit Windrauschen. Im Output wird das geschätzte Rauschen dann durch die Cepstrum-Koeffizienten von Referenzwerten ersetzt, um das Rauschen danach herausfiltern zu können.
Mein Code zur Bearbeitung des Signals sieht im Groben so aus:
Code:
% aktueller Abschnitt des Signals (inkl. Fensterung)
signal_in = signal(cnt:cnt+block_size-1).*hann_win;
signal_in_fft = fft(signal_in,fft_size);
signal_in_fft = signal_in_fft(1,1:bins); % redundanten Teil entfernen
Für mein Einsatzgebiet reicht bei einer fft_size = 512 ein threshold von 30. Der nichtredundate Teil der FFT ist bin = fft_size/2+1.
Wenn ich da nun einen threshold von 0 wähle, besteht das Spektrum signal_out_fft nur aus dem Eingangsspektrum und in der Filterung wird alles komplett entfernt (da Eingang und Ausgang gleich sind); das Signal scheint also richtig wiederhergestellt zu werden.
Einstellungen und Berechtigungen
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