guidance for designing a pilot signal Introduction
I wish to design a pilot signal which will be used to estimate the channel, and will be sent prior to data, on a wireless acoustic link (underwater acoustic communication).
I've been reading on the matter, but I'd like some pointers. I should state that I don't care about things like bandwidth regulatory restrictions and such as no restrictions exist. I'm however limited by the bandwidth of my analog frontend, and possible sampling frequency.
I intend to use a matched filter at thr receiver which hopefully will be able to provide me the delays between echos.
My thoughts so far
I thought of using chirps due to their natural resilience to multipath. I currently have the following questions:

I consider using a single longer chirp, or a few shorter ones, possible "coding" them by using a short pn sequence. Is there an obvious preference to use either one?
I considered using some pulses shaping filter, so as to to incorporate it into my match. filter. My fear is that without it the analog frontend will distort the signal affecting the matched filter performance. Am I correct to consider this? How will the pulse shaping affect the autocorrelation properties of the chirps.
If i use a root raised cosine filter, are there any standard guidelines regarding choice of beta?

answers as well as relevant links will be appreciated

An ideal pilot signal for channel estimation (without any apriori further knowledge of the channel) would transmit power uniformly across the channel and be easy to replicate in the receiver for channel estimation, as channel estimation cannot be determined where no signal exists for comparison, as well as carrier and timing recovery for a synchronous digital modulation data link. PN sequences are good choices due to their characteristics that mimic white noise. Further and due to the same properties resulting in an auto-correlation that appears as a single impulse as expected with white noise, they are also useful for synchronization such as estimating the start of a data packet. The PN sequence also provides fast convergence for timing and carrier recovery in a digital data link, but you may also consider adding a certain duration of a 10101010 sequence after the PN sequence (if the overhead allows) which will provide further robustness for timing recovery.
A chirp could also be used providing both properties of having a white spectrum (on average across the frequency band of interest for which the pulse is applied) and a sharp autocorrelation function (as demonstrated by the matched filter output being a “compressed pulse) although the use of PN codes could more easily offer a greater set of uncorrelated patterns that could be used in the same pilot time interval which could be used for other features related to packet identification (if desired). A chirp does not readily provide for carrier recovery since it is difficult to distinguish carrier offset from time offset, nor does it provide information on the data rate that the toggling of the PN sequence would provide.
The primary use and purpose of a pulse shaping filter is for applications where the spectrum must be contained. (With no pulse shaping filter, and the transmission of rectangular pulses, the ideal matched filter becomes and integrate and dump filter (or equivalent moving average process) and there is no penalty in doing this other than the much wider transmit spectrum. Distortion to pulse shaping is a non-linear distortion mechanism such as  driving amplifiers too hard into compression rather than linear distortions such as pass-band ripple, reflections etc that can be corrected for with channel equalization (as these distortions end up being similar to multipath that the channel introduces).   The resulting spectrum after a pulse shaping filter is also flatter across a narrower spectrum (just over B wide for a symbol rate of B symbols/sec) while a rectangular pulse will have about a 4 dB roll-off over that same bandwidth. For this reason the use of a pulse shaping filter for channel estimation patterns will result in more efficient channel estimation applying more equal power across the channel. When timing is of primary concern (such as with GPS), pulse shaping will degrade time resolution since for that the widest bandwidth is desired (for the sharpest autocorrelation peak). To create a matched filter, it is typical to decompose the pulse shaping into the cascade of two identical filters (thus square root filters), with one in the transmitter and the other in the receiver.
Here is a link that describes the trade space with pulse shaping further, but again first consider what your spectrum containment needs to be and then design from that; the tighter the pulse shaping the more complicated the receiver will be both with the duration of the matched filter and in the timing recovery algorithm/loop.
Pulse shaping with RRC : Number of taps


Definition of minimum-phase system I saw a couple of definitions for minimum-phase in different textbooks and I'm trying to understand what the implication of each of them.
The first definition I saw was: An invertible system which both it and it's inverse are causal and (BIBO) stable.
The second definition is: An invertible system which both it and it's inverse are causal and have finite energy.
It was written in one textbook that had the second definition as an alternative to the first one that these constraints allow the system to have zeros and poles on the unit circle and I don't understand why (as opposed to the first def).
just to make sure, I'm not restricting myself to rational systems.
Another thing I saw in some textbook and didn't understand is that the second definition implied that the transfer function of the system and it's inverse are analytic in the exterior of the unit circle.
So, to conclude, my questions are:

why we have both of the definitions of minimum-phase?
why the 2nd def implies that the system can have zeros and poles on the unit circle as opposed to the first def?
why the 2nd def implies that the system and it's inverse have analytic transfer functions in the exterior of the unit circle?

Sorry for any grammar mistakes and thanks for any clarifications.

one thing about a non-minimum phase system (with a rational transfer function), is that it can be thought of as the series concatenation (or cascade) of a minimum-phase system, having identical magnitude response as the given non-min-phase filter, with an all-pass filter.  the APF will have a poles that cancels specific zeros of the min-phase system that are inside the unit circle, and puts in new zeros ("reflected zeros") outside the unit circle.  so the new system always has more phase shift than the min-phase system.

so here's another definition of a minimum phase system:

An LTI system or filter

$$ H(f) = |H(f)| e^{j \arg\{H(f)\}} = |H(f)| e^{j \phi(f)} $$

is minimum phase if and only if the natural phase response, in radians, is the negative of the Hilbert transform of the natural logarithm of the magnitude response:

$$ \phi(f) \triangleq \arg\{ H(f) \} = -\mathscr{H}\big\{ \ln( |H(f)| ) \big\} $$

since

$$\begin{align}
 H(f) & = |H(f)| e^{j \phi(f)} \\
  & = e^{\ln(|H(f)|)} e^{j \phi(f)} \\
  & = e^{\ln(|H(f)|) + j \phi(f)} \\
  & = e^{\ln(H(f))} \\
\end{align}$$

this is relating the real and imaginary parts of the complex natural logarithm of the frequency response.

$$ \Im\{ \ln(H(f)) \} = -\mathscr{H}\big\{ \Re\{\ln(H(f)) \} \big\} $$

so a question to ask, why is the pole/zero definition stated in the question equivalent to this Hilbert transform definition??


Bandwidth of the Derivative of a Message If I have a message signal $m(t)$ and it has a bandwidth $B$. I know that the bandwidth of $m^N(t)$ is $NB$. But what is the bandwidth of $\frac{d m(t)}{dt}$?
Thanks!

Hint:

Checking a table of Fourier transform properties you'll find that the Fourier transform of the derivative is given by

$$\mathcal{F}\left\{\frac{dm(t)}{dt}\right\}=j\omega M(\omega)\tag{1}$$

where $M(\omega)$ is the Fourier transform of $m(t)$. Now if $M(\omega)$ is band-limited, what does this tell you about the Fourier transform of $dm(t)/dt$?


Are there alternatives to the bilinear transform? When designing a digital filter based on an analog filter we usually use the bilinear transform. To approximate a discrete transfer function $D_a(z)$ from analog (continuous) transfer function $A(s)$ we substitute

$$z = \frac{1+sT/2}{1-sT/2}$$

where $T$ is the sampling period. Alternatively, to approximate an continuous transfer function $A_a(s)$ from discrete transfer function $D(z)$ we substitute

$$s = \frac{2}{T} \frac{z-1}{z+1}$$

Are there alternative methods of performing such conversions? Are there better approximations?

Analog filters are stable if the poles are in the left half of the s-plane (figure on the left) and digital filters are stable if the poles are inside the unit circle (figure on the right). So mathematically all that is needed to convert from analog to digital is a mapping (conformal?) from the half-space to the unit disk and the $\jmath\Omega$ axis to the unit circle $\vert z\vert=1$. Any transformation that does this is a possible candidate for being an alternative to the bilateral transform.

Two of the well known methods are the impulse invariance method and the matched Z-transform method. Conceptually, both of these are similar to sampling a continuous waveform that we're familiar with. Denoting the inverse Laplace transform by $\mathcal{L}^{-1}$ and the Z transform as $\mathcal{Z}$, both these methods involve calculating the impulse response of the analog filter as
$$a(t)=\mathcal{L}^{-1}\{A(s)\}$$
and sampling $a(t)$ at a sampling interval $T$ that is high enough so as to avoid aliasing. The transfer function of the digital filter is then obtained from the sampled sequence $a[n]$ as
$$D_a(z)=\mathcal{Z}\{a[n]\}$$
However, there are key differences between the two.
Impulse invariance method:
In this method, you expand the analog transfer function as partial fractions (not in the matched Z transform as mentioned by Peter) as
$$A(s)=\sum_m \frac{C_m}{s-\alpha_m}$$
where $C_m$ is some constant and $\alpha_m$ are the poles. Mathematically, any transfer function with a numerator of lesser degree than the denominator can be expressed as a sum of partial fractions. Only low-pass filters satisfy this criterion (high-pass and bandpass/bandstop have at least the same degree), and hence impulse invariant method cannot be used to design other filters.
The reason why it fails is also quite clear. If you had a polynomial in the numerator of the same degree as in the denominator, you will have a free standing constant term, which upon inverse transforming, will give a delta function that cannot be sampled.
If you carry out the inverse Laplace and forward Z transforms, you'll see that the poles are transformed as $\alpha_m \to e^{\alpha_m T}$ which means that if your analog filter is stable, the digital filter will also be stable. Hence it preserves the stability of the filter.
Matched Z-transform
In this method, instead of splitting the impulse response as partial fractions, you do a simple transform of both the poles and the zeros in a similar manner (matched) as $\beta_m\to e^{\beta_m T}$ and $\alpha_m\to e^{\alpha_m T}$ (also stability preserving), giving
$$A(s)=\frac{\prod_m (s-\beta_m)}{\prod_n (s-\alpha_n)}\longrightarrow \frac{\prod_m \left(1-z^{-1}e^{\beta_m T}\right)}{\prod_n \left(1-z^{-1}e^{\alpha_n T}\right)}$$
You can easily see the limitation of both these methods. Impulse invariant is applicable only if your filter is low pass and matched z-transform method is applicable to bandstop and bandpass filters (and high pass up to the Nyquist frequency). They are also limited in practice by the sampling rate (after all, you can only go up to a certain point) and suffer from the effects of aliasing.
The bilinear transform is by far the most commonly used method in practice and the above two are rather more for academic interests. As for conversion back to analog, I'm sorry but I do not know and can't be of much help there as I hardly ever use analog filters.


How to reconstruct the signal if there were some temporary offsets I am not a signal processing expert. I have data like this:

where the red is what I have and the black is what I want to have. Any idea how to reconstruct this? The code for generating the data is here:

close all
Ndata = 1000;
X = nan(Ndata,1);
sigma = 0.5;
Y(1) = 50;
offset = 0;
for i =2:Ndata
    % consumption
    X(i,1) = -0.8*rand;
    if rand<0.01 % supply
        X(i,2) = 50;
    end
    if rand<0.5 % offset change
        if abs(offset)==0
           offset = exprnd(50)*(randi(2)*2-3);
        else
           offset = 0;
        end
    end
    X(i,3) = offset ;
    % noise
    X(i,4) = sigma*randn;
    Y(i) = Y(i-1)+sum(X(i,1:2));
    Z(i) = Y(i)+sum(X(i,3:4));
end
plot(Z,'r')
hold on
plot(Y)

The tricky thing is the fact that there is not only a small white noise, but there are also temporary excesses that are relatively large.

Thank you for any advice.

I just had a go at something stupidly simple: applying a median filter.

It's not ideal, but it might be a pointer to how to proceed.

The top plot is my instance of your plot.  The bottom plot shows the "noiseless" line from the top plot, with the median filtered version in black.

N2= 10;
for k=1:1000,
    idx1 = max(1,k-N2);
    idx2 = min(1000,k+N2)
    M(k) = median(Z(idx1:idx2));
end


Low pass fir filter with non-unity gain In my filter specs, the passband gain is given as 6.02db. How do I achieve this?

Usually all texts give the passband gain as 0db (as well as passband ripples if any).

As pointed out by niaren, you have to multiply by two, but this also means that your whole magnitude response is shifted up by 6dB, also in the stopbands. So you will have 6dB less stopband attenuation. Of course, the attenuation relative to the passband gain remains constant. You have to take this into account when designing your filter. Say you need 60dB stopband attenuation in absolute terms, not relative to the passband gain, then you need to design a 0dB gain filter with 66dB stopband gain. After multiplying the coefficients by 2 you'll get 6dB passband gain and 60dB stopband attenuation.


Setup frequency array properly I've just stepped into the FFT's world since the last weeks, therefore I still don't have an clear overall idea of how to setup the frequency array.

In case of even number N, I usually do (with respect to the Nyquist cut-off frequency)

$f=\left(-N/2\,:\, N/2-1\right)*f_{s}$

but how does it work in case of odd number N. Is it just

$f=\left(-(N+1)/2\,:\,(N+1)/2-1\right)*f_{s}$

I thank you all in advance.

Posted for anyone who may find this useful...

I created a picture that shows DFT frequency bin spacing for odd and even cases of N where N is the number of samples.  FFTs usually operate on an even number of samples (the algorithm works by repeatedly breaking the problem into halves), so only the even case applies.  The DC component (0*fs) is always part of the computation .  The bin at the nyquist frequency (fs/2) only shows up for the even case.  The frequencies above the nyquist frequency reflect back on the negative side of the frequency axis in reverse order.  For real signals the reflected FFT values are complex conjugates of their positive counterparts.  There are two specific cases depicted on the complex plane, N=5 for the odd case and N=8 for the even case.


Determine spectrum amplitudes for half-wave rectified sine I am trying to learn how to solve a bunch of digital signal problems and I have trouble understanding the solutions provided by this book I'm using. Basically, this problem asks me to determine spectrum amplitudes for half-wave rectified sine, as the title mentions. The signal looks like this:

The development in complex (exponential) Fourier series will be given by:

What properties are used to get to this answer? I don't get from where did the
 came from.

Full solution:

I was missing the red n. From that and using

solved it.

You have to integrate by parts twice. Leaving out constants you have

$$\int_0^{T/2}\sin(\omega_0t)e^{-jn\omega_0t}dt=-\frac{\cos(\omega_0t)}{\omega_0}e^{-jn\omega_0t}\bigg|_0^{T/2}-jn\int_0^{T/2}\cos(\omega_0t)e^{-jn\omega_0t}dt$$

If you again perform integration by parts on the integral on the right-hand side of the above equation you get the final result:

$$\int_0^{T/2}\sin(\omega_0t)e^{-jn\omega_0t}dt=\frac{1+(-1)^n}{\omega_0}+n^2\int_0^{T/2}\sin(\omega_0t)e^{-jn\omega_0t}dt$$

If $I$ is the integral on the left-hand side you have

$$I=\frac{1+(-1)^n}{\omega_0}+n^2I$$

and finally

$$I=\frac{1+(-1)^n}{\omega_0(1-n^2)}$$

Note that you have to consider the case $n=1$ separately. Also note that you could also compute the coefficients $c_n$ by direct integration using $\sin(\omega_0t)=(e^{j\omega_0t}-e^{-j\omega_0t})/2j$.


How to find frequency response, stability, and causality of a linear system? I have the following transfer function:

$$H(s)=\frac{s}{(s+1)(s+2)}$$

How can I find the gain and phase response of the above system? I know the first step has something to do with substituting $s = j\omega$ into $H(s)$. How can I find whether the system is stable and whether it is causal?

NOTE: This answer is community wiki: If you feel there is a need for
  correction -please update the answer right here.

1. Frequency response:
For any system - the (discrete time or continuous time) Fourier Transform of the impulse response is same as the frequency response of it. You can replace $s$ by $j\omega$ to work on this.

2. Stability:
System is stable if, the bounded input produces bounded output. The definition of can be found here in the wiki link.  For this, the necessary condition to prove is to see if the impulse response be absolutely integrable, i.e., its L^1 norm exist.

For a rational and continuous-time system, the condition for stability is that the region of convergence (ROC) of the Laplace transform includes the imaginary axis.

3. Causality:
Any system who's impulse response at sample $s_i$ doesn't require to know any samples presented $s_{i+1}$ or after wards, than the system is causal.

EDIT:
Now given the above i am trying to put the calculation as per your system.

$$ H(s)=\frac{s}{(s+1)(s+2)}  $$

Now we can replace $s$ with $j\omega$ to get the frequency response in terms of magnitude and phase response.

$$ H(j\omega) =  \frac{j\omega}{(j\omega+1)(j\omega+2)} $$
$$ H(j\omega) =  \frac{j\omega}{2 - \omega^2 + 3j\omega } $$
$$ Magnitude =  20 Log_{10} \sqrt { ( H(j\omega) )^2 } $$
$$ Magnitude =  20 Log_{10}  \sqrt { \frac { (j\omega)^2 }{(2 - \omega^2)^2 + (3j\omega)^2} }  $$

NOTE: square is applied separately on real terms than on imaginary terms! I got this understanding from this document here.

Replacing $-j^2$ by $-1$

$$ Magnitude =  20 Log_{10} \sqrt { \frac{-1\cdot\omega^2}{(2 - \omega^2)^2 - (3\omega)^2} } $$

$$ Magnitude =  20 Log_{10} \sqrt { \frac{-1\cdot\omega^2}{(4 - 5\omega^2 +\omega^4)} } $$

$$ Phase = arc tan (real / imaginary)  $$

$$ Phase = -1 \cdot \tan^{-1}(\frac {-3\omega}{(2-\omega^2)} )   \ $$

Please correct me if i am wrong.

Another literature that can help you.


Getting the amplitude and phase of harmonics in C I have a periodic series of 360 values. How can I use the fftw c libraries to get the amplitude and phase of the first, second, third and fourth harmonic?

If I do this for N=360

fftw_plan_r2r_1d(N,input_array,output_array,FFTW_R2HC,FFTW_ESTIMATE);

I can get the 360-point Discrete Fourier Transform (DFT) where output_array is complex Hermitian in half-complex format.

Where k$^{th}$ component of output_array is the real part and (N-k)$^{th}$ component of output_array is the imaginary part.

So if I want the amplitude or magnitude(??) of the first harmonic should I do this?

ampl_1sth=sqrt(output_array[1]*output_array[1]+output_array[N-1]*output_array[N-1])

Is that correct? then the amplitude of the second harmonic will be the same but with output_array[2] and output_array[N-2] and so on right?

And then how can I get the phase of the first, second, third and fourth harmonic with the info from the DFT?

For a reference to fftw see here www.fftw.org/fftw3.pdf

Thank you.

You just need to know how to compute magnitude and phase of a complex number $z=z_R+iz_I$:

$$|z|=\sqrt{z_R^2+z_I^2}\\
\arg\{z\}=\left\{\begin{array}{ll}0,&z_R=z_I=0\\
\pi/2&z_R=0,\; z_I>0\\
-\pi/2&z_R=0,\; z_I<0\\
\arctan\left(\frac{z_I}{z_R}\right),&z_R>0\\
\arctan\left(\frac{z_I}{z_R}\right)+\pi,&z_R<0,\;z_I\ge0\\
\arctan\left(\frac{z_I}{z_R}\right)-\pi,&z_R<0,\;z_I<0\\
\end{array}\right.
$$

Luckily, the computation of the phase $\arg\{z\}$ is simplified by the existence of the function atan2 in most programming languages. Watch out, you need to use it like this:

$$\text{phase = }\text{atan2}(z_I,z_R)$$


Find the Fourier transform of $g(k)$ from $G(z)$ for frequency $=1/2$ $$G(z)=\displaystyle \frac{\frac{1}{z}}{1+\frac{5}{6z}+\frac{1}{6}z^{-2}}$$

I found:
$$g(k)=\displaystyle \left(\frac{-1}{3}\right)^k - \left(\frac{-1}{2}\right)^k$$

I don't understand how I can find the fourier transform of $g(k)$ from $G(z)$.

First of all it's very important to realize that without any further information it is not possible to uniquely identify the corresponding time-domain sequence from the given $G(z)$. The reason is that you need to know the associated region of convergence (ROC). If the ROC is not known then for the given $G(z)$ there are three possible time-domain sequences. One of them is causal, another one is anti-causal, and the third one is two-sided.

If the ROC is given as $|z|>\frac12$, or, equivalently, if you know that you're looking for a causal sequence, then the inverse transform is uniquely defined. It is similar to the one you found, but there is an important difference:

$$g[k]=\left[\left(-\frac13\right)^k+\left(-\frac12\right)^k\right]u[k]$$

You need the unit step $u[k]$, which makes the sequence causal. The sequence you found extends to infinity for positive as well as for  negative indices $k$, and does not even have a $\mathcal{Z}$-transform.

Now that you know that the ROC includes the unit circle, you know from the answers to this question that in order to find the Fourier transform you can simply evaluate the $\mathcal{Z}$-transform on the unit circle, i.e. set $z=e^{j\omega}$. What they didn't tell you is that this is only possible because the unit circle is included in the ROC. Otherwise you can't do that because the Fourier transform wouldn't even exist.


FFT and Power Spectrum Normalization On many websites, including MathWorks, it was suggested to normalize the fft spectrum (MATLAB or numpy) by dividing it by the total number of samples ($N$). For a sinusoidal signal, for example:

$$x(t)=5cos(2\pi f_0t)$$

This produces a two-sided spectrum peak at $f_0$ with a peak amplitude of 2.5. However, that's not the result of the exact fourier transform, which should be:

$$X(f)=2.5(\delta(f-f_c)+\delta(f+f_c))$$

That is, I expect the amplitude to increase as the frequency bin shrinks (with the total area remaining at 2.5). That is not the case if I normalize with respect to $N$. If I normalize instead with respect to the sampling frequency ($F_s$), then it works out for the sinusoidal signal. The total power density is obviously finite ($P=\frac{1}{T}\int_{0}^{T}x(t)^2dt=12.5$); but the numerical power spectrum (unit^2/Hz) will have varying results depending on the bin size.

Questions:

Can I assume that the division by sampling frequency normalization would be more physically relevant for vibration analysis, if I want to extract things like $g^2/Hz$?
If my signal is a finite sum of sines, what's the best way to visualize the power? I would expect the power spectrum to vary again with the total sample time (i.e. not very useful).
If I assume all continuous signals can be represented by fourier series, how come non-sinusoidal signals have finite and well-defined power spectrum?

I would say most of your confusion arises from conflating the continuous case with the discrete case of the FT.  The are similar, and related, but not equivalent.  "Area" is a continuous concept.  "Sum" is discrete.

I normalize my DFTs by 1/N principally because it makes the bin values essentially sample count independent.  Other reasons too.  For the sum of squares to hold, you should normalize with $1/\sqrt{N}$.

It would be helpful to know if you have a Linear Algebra background in order to better explain some of these things.  Also if you are talking about the continuous case or the discrete one.

How much energy is there in the fluctuation of a color pattern in a fabric?  Meaningless question, right?  The concept of energy is an application issue and is not inherent in the DFT.  What is inherent is the ability to sum squares of values.  With a $1/\sqrt{N}$ normalization, the discrete fourier transform can be represented as the multiplication of a unitary matrix, thus the sum of squares is preserved.

$$ Z = F S $$

Suppose Z is a complex vector of the DFT bins, F is the tranformation matrix, and S is a complex vector with your signal.  Of course, a real valued signal will fit too.

Let ($^*$) stand for Hermitian, aka, the conjugate transpose.

$$ Z^* = S^* F^* $$

$$ Z^* \cdot Z = ( S^* F^*  ) \cdot F S $$

$F^*$ also happens to be the inverse DFT so $ F^* \cdot F = I $, the identity matrix.

$$ Z^* \cdot Z = S^* \cdot I \cdot S  = S^* \cdot S $$

The left and right parts of these equations are simply the sum of the squares of the magnitudes.

Between the forward and inverse DFTs, a factor of $1/N$ needs to be included.  To split it evenly, it becomes $1/\sqrt{N}$ in each direction.

The "energy" term is associated with the sum of squares.  The units are dependent on the application.  The DFT doesn't care.


Fixed Point Design Resources Can anyone point me to good methodologies for designing fixed point versions of possibly nonlinear signal processing algorithms?
Are there any systematic methods other than simulation for optimizing and analyzing finite word length effects?  References would be appreciated!

Two algorithms of interest are IQ balancers and digital AGC loops.

Make your block diagram
At each point where significant quantization can happen, add noise
Analyze your system's behavior with that added noise

If you know that the quantization effects will be essentially random, and if the following stages tend to low-pass or band-pass filter, then model the quantization noise as Gaussian with $x_n \sim N(0, q/12)$, where $q$ is one LSB.

The justification for this is that the quantization noise will be uniformly distributed (hence $\sigma = q/12$), and the following filter's output will contain a weighted sum of a bunch of samples of the quantization noise, which will tend to Gaussian by the central limit theorem.

If you want to be Maximally Paranoid, then figure out what the worst-case behavior of the quantization noise would be for your application (stuck high, stuck low, some specific signal, etc.), and model the quantization noise as having that shape, with a magnitude of $q$.  In the case of a linear system that has any frequency selectivity, this works out to a square wave at a frequency equal to the systems highest sensitivity to noise at the injection point (or a sine wave with an amplitude of $1.09 q/2$.

The justification for this is really just paranoia -- but if you have an absolute upper bound on the effects of quantization, and the system still works well enough, then you just know you're done.

To design data path widths, you can turn this around, and determine the acceptable level of quantization at each step, and make sure that your data paths are wide enough (and scaled correctly) so that quantization is smaller than your acceptable level.

It's dangerous to make general statements about nonlinear systems, but I would hazard a guess that if you can analyze the algorithm on paper at all, there's a good chance that you can analyze the algorithm plus quantization on paper.


Band Limited Impulse Train Synthesis I'm currently working on a piece of DSP/Synthesis software and I have run into the issue of generating band limited waveforms. I have been doing some research to try and find a good solution but I haven't had much luck. I have found some things about a BLIT to generate complex waveforms but I have not found a good example of an implementation. I am not very good with Math or Mathematical notation but I am great with dealing with c/c++ code. If anyone has a good explanation of how to generate a BLIT using c/c++ or resources about that I would much appreciate any help that you could give.

Thanks

P.S. I know this could also be a Stack Overflow question because of the programming nature but I thought I would try the DSP community first. Thanks again.

the main idea behind BLIT is that these analog synth waveforms that we are trying to generate digitally can be thought of as the integral (over $t$) of impulse trains.

a sawtooth can be thought of as the integral of the sum of a little bit of DC and an impulse train.  a square wave is the integral of impulses of alternating signs.  the  triangle wave is the integral of the square wave.

so, to create bandlimited waveforms of the above, the impulse trains are bandlimited which means that each impulse $\delta(t-t_n)$ is replaced by a $\operatorname{sinc}(t-t_n)$ function, which is that impulse bandlimited through a Nyquist brick-wall LPF.  that sequence of bandlimited impulses is a BLIT.

then, since integration is a filter with s-plane transfer function of $H(s)=\frac{1}{s}$ and is LTI (Linear, Time-Invariant), integrating the BLITs will introduce no new frequency components.  if your BLITs are bandlimited, so are the other waveforms that are derived from filtering the BLITs.


What Does Multi Scale Image Analysis Mean? I'm reading an article about multi-scale representation of image, and it is said that, convolving the image with a Gaussian kernel at different $\sigma$, then different scale representation is created.

This really confuses me, because I think the image is just blurred by Gaussian with different $\sigma$, each blurred version still has the same image size, so I don't think they're scaled, unless each blurred is downsampled to a smaller image size, which means I think different size representation is truely the

Multi Scale means different scale of details of the image.
It comes from the intuition that blurring image with Gaussian Blur with $ \sigma $ vanishes data which is smaller then $ \frac{\sigma}{3} $.

In practice it is similar to Multi Rate Analysis of Filters in Signal Processing.
You look at the image in frequency domain and divide its data into octaves.
This partition can be done with any LPF Filter (Yet for perfect reconstruction there is some analysis).
In Image Processing the most widely used LPF is the Gaussian Blur.

By the way this is the idea behind Wavelets as well (Intuitively, in Wavelets they create special family of filters with certain set of properties).


Timing synchronization in limited channel I need measure timing of an signal that is send through a noise band limited channel (wireless).
By measure timing, I mean measure as exactly as possible when the signal has arrived at the receiver.
My idea is to use a matched filter (correlator) to match some symbol send by the transmitter and then locate the maximum point to find the exact time of arrival.

But is this the best way to do it?

There are numerous approaches to make timing synchronization at the receiver. Yes, they are based on correlation in general.

But some rought methods also exist, e.g. you can track power of the input signal by two windows, and comparing measured window power values, you can detect begining of the signal. Of course it will work unreliably at bad SNR scenario. And also some time uncertainty will exist.

So the matched filter is quite good idea. If you are dealing with streaming modem, you can use partly signal periodicity to detect signal at first and then make timing (like it's done in DVB-T for example). If packet modem is your case, you have to use frame header for the detection and timing problems.

The main idea here that if you work with stream you can loose some data, while you are in synchronization process (you can trade off synchronization rapidness and its hardware complexity). But if you work with packet communications you're not allowed to loose frames if it's possible. So you must complete your synchronization with the aid of frame header only.

And if you're not sure about the type of sync symbol yet, try maximum length sequences for example, they are quite common for such a task.

Hope this helps.


.ptw Video conversion On Video recording with IR camera, a file in .ptw is obtained. I want to analyze the file using standard software, like MATLAB. However, video processing in MATLAB requires the file to be converted to some standard formats like .avi or .mp4 or .wmv files. Please provide a way to make similar conversions, with the help of a software tool or any matlab function etc.

Unfortunately the only way I know of to convert their .ptw format to an standard video format is by using FLIR's ResearchIR software:  http://www.flir.com/thermography/americas/us/view/?id=51371

ResearchIR also supports exporting to .mat files too.

Hope that helps,
Andy


Frequency range for STFT in Librosa I am currently looking at the STFT for Librosa. I was wondering how to understand the output of the STFT function, specifically what kind of frequencies the different values represent.

Say I have a n_fft of 256, that means that the shape output of the STFT will be 129 (1+n_fft/2). I therefore understand that for each frame I have 129 bins, these bins represent some sort of frequency within this frame. I am assuming that the bins start with a frequency of 0hz and go up to some value, with each bin representing an equally large range.

But how do I figure out what the maximum frequency used for STFT is? And what about the range of each of the 129 bins?

Just to be clear, I have tried looking at the source-code and at the documentation, but I have not really become much wiser.

Having a look into the sources of librosa.display.specshow reveals how bins are converted into frequencies internally: The plotting uses librosa.core.fft_frequencies, which shows that it is basically the same as numpy.fft.fftfreq:
>>> librosa.fft_frequencies(sr=22050, n_fft=16)
array([     0.   ,   1378.125,   2756.25 ,   4134.375,
         5512.5  ,   6890.625,   8268.75 ,   9646.875,  11025.   ])


Sampling at subnyqvist frequencies I went through this material on Sampling theorem. http://www.wescottdesign.com/articles/Sampling/sampling.html

It uses a powerline signal example and concludes that if a signal is cyclic and repeating over a period N then you could sample the signal at lesser than N samples per second.In the example he shows for a 60Hz power signal he can sample and recover the signal by sampling at 19 Hz. He uses an expression

Ts = (M + P/N )*1/F where M>0 and N>0 where p is not equal to 0. I haven't understood it. Can anybody here tell me how it is done.

The point in the article is that you can sample the signal at lower frequencies than 60Hz by subsampling. Consider the following two scenarios of sampling the 60Hz signal:

You sample like Nyquist wanted, at freq $\geq 120\mbox{Hz}$. In one period of your 60Hz signal, or $T=16.67\mbox{ms}$, you will have a correct representation of your signal via sampling.

You sample at a frequency even lower than the $60\mbox{Hz}$ signal freq. The trick is to sample at frequencies that do not divide the main frequency, which means you cannot use $2, 3, 6, 10, 12, 15, 20, \mbox{ or } 30 \mbox{Hz}$ frequencies, because then you'd always get the same point of the signal.

Let's say you sample at $19\mbox{Hz}$ like in the article. The whole point of this is to save the points on your graph, and get the desired signal point-by-point. So you start sampling your signal, get one point, then whole 3 periods of the signal will pass, and you get a second point, which you put on the graph, and so on until you have enough.

How many points do we need? The Nyquist theorem you stated tells us the number of points we need to be sure we got 100% of the signal. $$T_s = (M + \frac{P}{N})\frac{1}{F},$$  where $F$ is our signal frequency of $60\mbox{Hz}$, $T_s = 1/F_s = 1/19\mbox{s} = 52.63\mbox{ms}$ is our sampling period, and $M, N, P$ are integers where $N$ and $P$ must have no common factors. If we try to get such integers, we can see that for $N = 19, M = 3, \mbox{ and } P = 3$ our equation holds. So, we need 19 points for our graph, in order to have the whole signal. This means we need to wait 19 times our sampling period to get our total time of sampling, so we get $$ T_{tot} = 19\cdot 52.63\mbox{ms} = 999.97\mbox{ms} $$

So we see that sampling is possible at even lower frequencies, but at much longer time :) Hope this clarified your question! :)


Nyquist Rate (Sampling Frequency) for $ {f}^{2} \left( x, y \right) $ We are given that $f(x,y)$ is highest frequency is $\omega$ what will be the frequency sample rate if we want to restore the function of the form $g(x,y)=f^2(x,y)$

Would it be correct to say that because $\sin^2x=1-\cos 2x$, we will have to sample it at $4\pi$ ?

By the Convolution Theorem multiplication in Time / Spatial domain is equivalent of Convolution in the Frequency Domain.

The sampling rate (In its classic interpretation) is proportional to the support of a function in the frequency domain.

So if a function has a certain support in frequency, what would be its support after convolution with itself? Indeed it will be doubled in each dimension.


what are resources you recommend to learn GPS signal processing? I want to learn about GPS and other GNSS signal processing algorithms.
I have solid background in signal processing and some experience in OFDM.
Could you please recommend some books or online resource?

A good place to start is Van Diggelen and Enge's old course, now on YouTube.
There's a lot of GNSS related information on Wikipedia but I would look at Navipedia first.
Finally, you could try one of open source software receivers, the most comprehensive being gnss-sdr.


Books on delay / lag compensation Which books do you recommend that provide comprehensive overview of delay compensation techniques in control systems?

Perhaps the standard is Ogata's "Modern Control Engineering".

I also recommend Franklin/Powell/Emami-Naeini's "Feedback Control of Dynamic Systems".

For a more basic introduction you might consider Nise's "Control Systems Engineering", but I'm not a fan of his writing style.


Optimizing Two Parallel IIR Filter Sweeps Got an audio related question here; So I have two IIR low pass filters that work in parallel to do a filter sweep. Basically:

filtered_sample = filterSweepA(xn) + filterSweepB(xn);

End result: it sounds beautiful! However, now I want to optimize it. I was thinking about combining the two IIR filters into one so I'd only have to call filterSweep only once but I'm pretty sure you can't do that because of the recursive properties...

I think I understand that FIR filters would be the way to go but I've searched the internets and I've seen nothing about FIR filter sweeps.

Is this possible/just difficult to do?
If this can be implemented how would I be doing this?

I don't really know how to figure out FIR filter coefficients so I guess that would be the place to start...

FYI I am controlling the start, stop filter frequency cutoffs and the sweep time (usually around 1.5 seconds).

Two (linear, time invariant) filters in parallel, $H_1(z)$ and $H_2(z)$, are equivalent to a single filter $H(z) = H_1(z) + H_2(z)$ (i.e. the sum of both transfer functions).

You don't give any information on how those filters are represented. Normally, they are represented in terms of the recursive ($a_k$) and non-recursive ($b_k$) coefficients (see Wikipedia for reference).

Now, if you are going to add both transfer functions, then you need to find the common denominator (which is the product of two polynomials, one from each filter), and more cumbersome yet is to find the new numerator polynomial coefficients.

If you are doing this automatically (since you adjust your filters periodically), I suggest you keep your algorithm as is.

Executing an IIR filter is not very CPU intensive since the order of the filters is normally not greater than 10 or so.

Note: if the filters were FIR, then you simply add the filter coefficients (there are no recursive coefficients (and no division of polynomials) to worry about). But FIR filters are normally much more CPU intensive than comparable IIR.


Spectrum of OFDM with raised cosine window I'm having trouble implementing OFDM with a raised cosine (RC) window in Matlab. I know how to generate an OFDM signal and how to show its spectrum, I just don't know how to generate the window extensions. Each OFDM symbol is extended by $TW$ samples at both ends to smooth the transitions between successive symbols, this is done mainly to improve the out of band spectrum and reduce the interference to adjacent channels. I'm hoping someone here knows how to do this.

This spectral shaping technique is applied in time domain, after adding the guard interval (GI). But I find it easier to add GI and the cyclic extension for windowing in one step. Let $x(n)$ be an $N$ subcarriers OFDM symbol without guard interval. Then $W + G$ samples are copied to the beginning accounting for guard interval and windowing samples. Additionaly, $W$ samples are copied to the end, also for windowing:
$$
y(n) =
\begin{cases}
x(n+N) & \text{for} & -G-W \leq n \leq -1 \\
x(n) & \text{for} & 0 \leq n \leq N -1 \\
x(n-N) & \text{for} & N \leq n \leq N + W -1
\end{cases}
$$
In the next step, the raised cosine function is applied to the first and last $W$ samples of $y(n)$, respectively. The windowing function $w(n)$ is given by:
$$
w(n)=
\begin{cases}
\cos^2\left( \frac{n+G+1}{W-1} \frac{\pi}{2}\right) & \text{for} & -G-W \leq n \leq -G-1 \\
1 & \text{for} & -G \leq n \leq N - 1 \\
\cos^2\left( \frac{n-N}{W-1} \frac{\pi}{2}\right) & \text{for} & N \leq n \leq N+W-1 \\
0 & \text{otherwise}
\end{cases}
$$
$w(n)$ is similar but not equal to the transfer function of a raised cosine filter often used as impulse shaper in single carrier transmission systems. The two differences are: (1) the raised cosine function is applied in time domain for OFDM systems and in frequency domain for single carrier systems and (2) the "flat top" is usually much longer for OFDM systems, whereas its length is in a fixed relation with the flanks' length, given by the roll-off factor, for single carrier systems.

Finally, the OFDM symbol including GI and spectral shaping is calculated by
$$
z(n) = w(n)y(n)
$$

When transmitting several OFDM symbols $z_i(n)$, two consecutive symbols overlap at $W$ samples. The discrete transmit signal $u(n)$ is therefore given by
$$
u(n)=\sum_{i=-\infty}^\infty z_i(n-i(N+G+W))
$$

The implementation in Matlab should now be straightforward by substituting $n$ with $n' = n + G+ W+ 1$ in the above equations. Leave a comment if not.


I have some questions about background subtraction with Gaussian Mixture Model (GMM) I'm learning about Gaussian Mixture Model (GMM) for background subtraction.. And have some questions about it..

In Chris Stauffer paper "Adaptive background mixture models for real-time tracking" Gaussian probability density function is
$$\mathcal{N}(X_t,\mu,\Sigma)=\dfrac{1}{(2\pi)^\frac{n}{2} \left|{\Sigma}\right|^\frac{1}{2}}\exp\left(-\dfrac{1}{2}(X_t-\mu_t)^T\Sigma^{-1}(X_t-\mu_t)\right)$$
where $n$ in multi dimension, is $n$ in 3-D which is pixel position in x,y, and intensity of pixel?
In the same paper, they said:
$$\Sigma_{k,t}=\sigma_k^2.I$$
Is $I$ is Indentity matrix?
They also said red,green, and blue pixel values are independent and have the same variances, if it independent why R,G, and B have same variances? shouldn't it have different variances because the value of each channel is different..
If I set the algorithm with number of distribution is 3, so it's only 1 distribution that chosen as background model and the rest is foreground model? Because they said the first B distribution are chosen as the background model, where
$$B=\arg \min_b(\sum\limits_{k=1}^b w_k>T)$$

Yes, the dimension of $X$ should be 3, one for each of the red, green, and blue channels. You can see this from their Figure 2, where they show plots of the distributions of red and green channels only (to motivate why they want to use a Gaussian Mixture Model in the first place --- the distribute of pixel values is multi-modal).
Yes, $I$ is the identity matrix.
Random variables can have the same variances and still be independent. The variance is part of the model. In this case what the model is saying is: the red, green, and blue pixel values vary without regard to each other (they are independent) but the amount that they vary is about the same as each other (they have the same variances).
Whether only 1 distribution is chosen when $K=3$ will depend on the threshold $T$, surely?  Yes, it is possible that only one distribution will be chosen.


Speech corpora for VAD I'm doing a student research project connected with VAD.
And I want to find speech corpora with added noise.
I found many references to the NIST Speaker Recognition Corpora. But unfortunately, all links on main site are already dead.

Maybe someone has this Corpora? Or someone can recommend a good Corpora with noise?

Look into the MUSAN corpus from the Kaldi people (Dan Povey, et al):

Corpus audio/text
Check the arxiv paper

  MUSAN is a new corpus of music, speech, and noise. This dataset is
  suitable for training models for voice activity detection (VAD) and
  music/speech discrimination. Our corpus is released under a flexible
  Creative Commons license. The dataset consists of music from several
  genres, speech from twelve languages, and a wide assortment of
  technical and non-technical noises. We demonstrate use of this corpus
  for music/speech discrimination on Broadcast news and VAD for speaker
  identification.


Band energy and Parseval theorem In Python, I consider a dummy signal as follows:

fs = 512
sig = np.random.standard_normal((fs,))

sig corresponds to 1s of data at 512Hz.

To check Parseval theorem on the whole signal, I do:

fft = np.fft.fft(sig)
f = np.fft.fftfreqs(sig.size, d=1. / fs)
m = np.sum((np.abs(fft) ** 2) / sig.size)
e = np.sum(sig ** 2)
print('diff = %1.5f' % (m - e))

This actually returns 0.0. Now, to compute the energy of the signal in a given frequency band (say 5Hz - 15Hz), I would filter the signal using a bandpass filter and compute its energy as above:

def _butter(lowcut, highcut, fs, order=5, btype='bandpass'):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype=btype)
    return b, a

def _filter(data, lowcut, highcut, fs, order=5, btype='bandpass'):
    b, a = _butter(lowcut, highcut, fs, order=order, btype=btype)
    y = filtfilt(b, a, data)
    return y

filtered_sig = _filter(sig, 5., 15., fs)
e_band = np.sum(filtered_sig ** 2)

My question is: can I use Parseval theorem to compute the energy of the signal (in the 5Hz-15Hz) band using fft as above? I naively tried:

mask = np.logical_and(np.abs(f) >= 5, np.abs(f) <= 15)
m_band = np.sum((np.abs(fft[mask]) ** 2) / sig.size)

But the quantities m_band and e_band are not equal.
Where did I make a mistake?

The band-pass filter is not ideal. You will still have some frequency content outside the 5 Hz-15 Hz band. For a sufficiently high order filter, the difference between m_band and e_band should be small enough.


Manually Windowing an IIR Filter in Matlab I am currently trying to window an analog IIR filter in order to obtain an FIR filter. Though I realize the fir1 function exists within Matlab, I was wondering to see if it could be done step-by-step to produce a notch filter at 0 Hz to eliminate any DC gains in the system. However, I quickly realized that I may not know what I'm doing :)

Firstly, I created a butterworth filter containing the characteristics that I'm looking for just for a comparison. Once the Chebyshev II analog filter was created, the magnitude and phase response was viewed. It doesn't show a bandstop at 0 Hz, but my guess is because it hasn't yet been windowed to show this result? I also wasn't too sure how to combine the window and analog filter, and would love it if someone pointed me in the right direction :(

Fs = 250;
n = 5;
Rs = 40;
Wn = [0.1*2/Fs, 0.2*2/Fs];
[b0,a0] = butter(3,Wn,'stop');
[b1,a1] = cheby2(n,Rs,Wn,'stop','s');
freqz(b0, a0, 4096, Fs);
w = triang(250);
y = filter(b1,a1,w);

In the classical world of digital filters, an infinite impulse response (IIR) filter is used for its computational efficiency while achieving supreme performance compared to a finite impulse response (FIR) filter.

Usually, implementing an IIR filter happens in the time domain, with an architecture of direct form-I, II or else which provides a solution of the linear constant coefficient difference equation (LCCDE) describing the filter.

$$
\sum_{k=0}^{N} a[k]y[n-k] = \sum_{k=0}^{M} b[k]x[n-k]
$$

which is typically called an $N^{th}$ order LCCDE.

However, when there happens a need for implementing the IIR filter using time domain convolution instead;

$$
y[n] = h[n] \star x[n]
$$

then this requires the practical availability of the impulse response $h[n]$ of the IIR filter which is infinite length by definition.

To obtain a finite length, FIR, approximation $h[n]$ to the IIR digital filter, first choose your analog prototype filter, then apply the necessary transforms to get the digital IIR equivalent filter coefficients $a[k]$ and $b[k]$, and compute the truncated impulse response $h[n]$ of this filter, from $a[k]$ and $b[k]$ to a length that's adequate for the purpose of application and finally apply a further windowing to $h[n]$ if that's necessary...

The following MATLAB/OCTAVE command will give you a length $N$ truncated impulse response from IIR coefficients $a[k]$ and $b[k]$ with zero initial conditions.

h = filter(b,a,[1 zeros(1,N-1)]);


DFT derivative property? Does it have one? The continuous variant does, $f'(t) \rightarrow j \omega F(\omega)$ - but $jkX[k]$ definitely isn't it for DFT.
To find it there must be a useful simplification of $\text{DFT}(x[n] - x[n-1])$, where $x[n] = \sum_{k=0}^{N-1} X[k] e^{j2\pi kn/N}$, while also accounting for the one dropped sample (x[0]) - and I cannot see such a simplification.

Note: $x[n] = ...$ above is for the definition of IDFT I use, rather than where the signal derives from. The signal should be assumed the most general possible; not periodic, can be complex or real, but must be finite.
Note2: see my answer for some important caveats to (and direct contradicting of) the accepted answer, and optionally the comments discussions below answers. Thanks to @CedronDawg for the detailed discussion.
Note3: "derivative" = finite difference (for sake of this question).

To get a detailed answer along the lines of what you propose, we need to be careful about the normalization used in discrete Fourier transform (DFT) and inverse discrete Fourier transform (IDFT):
$$\text{DFT: }\quad X[k] = \sum_{n=0}^{N-1} x[n] e^{-j 2\pi kn/N}\tag{1}$$
$$\text{IDFT: }\quad x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j 2 \pi k n / N}\tag{2}$$
Those have a normalization that is directly compatible with fft and ifft from MATLAB, Octave, NumPy and SciPy. The indexes $k$ and $n$ run from $0$ to $N-1$. (MATLAB and Octave have a different indexing convention, $1$ to $N$.) Then:
$$\begin{array}{c}y[n] = x[n] - x[\operatorname{mod}(n-1, N)]\\
\begin{align}\\
Y &= \operatorname{DFT}\big(y\big)\\
&= \operatorname{DFT}(x*[1, -1, 0, 0, \ldots])\\
&= \operatorname{DFT}(x)\times\operatorname{DFT}([1, -1, 0, 0, \ldots])\\
&= X\times\operatorname{DFT}([1, -1, 0, 0, \ldots])\end{align}\end{array}\tag{3}$$
$$\Rightarrow Y[k]= X[k]\left(1 - e^{-j2\pi k/N}\right),\tag{4}$$
where $\operatorname{mod}$ gives the unsigned remainder, for example $\operatorname{mod}(-1, N) = N - 1$, the symbol $*$ denotes length-$N$ circular convolution and $\times$ denotes multiplication, and all sequences are of length $N$.
Circular convolution in the discrete time domain is exactly equivalent to multiplication in the discrete frequency domain, when DFT and IDFT are used to transform the sequences between the domains. See for example Circular Convolution - MIT OpenCourseWare. For "linear" convolution of discrete sequences, there is no such elegant pair of equivalent operations, which makes me think no expression proposed as an answer to the question will ever be fully satisfactory.
Considering input signals that extend to the left of the output range $0 \le n < N$ of IDFT, calculating the backward difference using modulo indexing is conditionally equal to calculating it without it:
$$x[n] - x[\operatorname{mod}(n-1, N)] = x[n] - x[n-1]\quad\text{conditionally}\tag{5}$$
under the condition that you only calculate it for some of the indexes:
$$0 < n < N,\tag{6}$$
or for $0 \le n < N$ under the condition that there's a hint of periodicity in the signal:
$$x[-1] = x[N-1].\tag{7}$$
A sufficient but not necessary condition that satisfies Eq. 7 is that $x$ is $N$-periodic, which is however prohibited by the signal defined as finite in the question. An example of another condition that satisfies Eq. 7 is $x[-1] = x[N-1] = 0$.
@CedronDawg's first answer provides a formula which corrects in the frequency domain the error in Eq. 5 if neither condition is satisfied: $Y[k] = X[k]\left( 1 - e^{-j2\pi k/N} \right) -  x[-1] + x[N-1]$. As an alternative derivation, in length-$N$ time domain the correction is an impulse:
$$\begin{align}&\big[x[0] - x[−1] - \big(x[0] - x[N-1]\big),\, 0,\, 0,\, \ldots\big]\\
= &\big[x[N-1] - x[-1],\, 0,\, 0,\, \ldots\big],\end{align}\tag{8}$$
which in frequency domain is a constant (see DFT Pairs and Properties: pair row 2, property linearity):
$$x[N-1] - x[-1],\tag{9}$$
to be added to all elements of $Y$ calculated by Eq. 4.
For a general $x$, the condition of Eq. 6 for Eq. 5 enables to use a length $N+1$ DFT and IFT to calculate the backward difference, by shifting the input to the DFT one step to the right, and finally by shifting the output from IDFT one step to the left. With forward difference $x[n + 1] - x[n]$, the shift would not be necessary, and I think this matches your discarding of the 0th sample. For a circular convolution implementation of convolution by a finite sequence, using a longer transform is a common trick to avoid the circular effects in a sufficiently large part of the output of the IDFT. For then obtaining the DFT of a partial IDFT output, in particular a one shorter, I don't think there is any shortcut.
We could perhaps re-express the question as: What is the DFT of the length $N-1$ forward difference of $x$ of length $N$, given $x$ and $X_{N-1} = \operatorname{DFT}(x_{N-1})$, a length $N-1$ DFT of the partial sequence $x_{N-1} = \big[x[0], x[1], \ldots, x[N-2]\big]$? Analogously to Eq. 4 we have:
$$y_{N-1} = x_{N-1}*[-1, 0, 0, \ldots, 0, 0, 1]\tag{10.1}$$
$$\Leftrightarrow Y_{N-1} = X_{N-1}\times\operatorname{DFT}([-1, 0, 0, \ldots, 0, 0, 1])\tag{10.2}$$
$$\Rightarrow Y_{N-1}[k] = X_{N-1}[k]\big(e^{j2\pi k / (N - 1)} - 1\big),\tag{10.3}$$
where each sequence is of length $N-1$. The desired forward difference $f_{N-1}$ is:
$$f_{N-1} = \big[x[1] - x[0],\, x[2] - x[1],\, \ldots,\, x[N-1] - x[N-2]\big].\tag{11}$$
Eq. 10.1 can be expanded to:
$$y_{N-1} = \big[x[1] - x[0],\, x[2] - x[1],\, \ldots,\, x[0] - x[N-2]\big].\tag{12}$$
By comparing Eqs. 11 and 12, it can be seen that:
$$f_{N-1} = y_{N-1} + \big[\ldots,\, 0,\, 0,\, x[N-1] - x[0]\big],\tag{13}$$
where the sequence in brackets is of length $N-1$. Taking the DFT of both sides of Eq. 13 and applying Eq. 10.3 gives the answer:
$$\begin{align}F_{N-1}[k] &= Y_{N-1}[k] + (x[N-1] - x[0])e^{j2\pi k/(N-1)}\\
&=X_{N-1}[k]\big(e^{j2\pi k / (N - 1)} - 1\big) + (x[N-1] - x[0])e^{j2\pi k/(N-1)}.\end{align}\tag{14}$$
This is the length $N-1$ DFT of the length $N-1$ forward difference of $x$ of length $N$.
Alternatively, you might be interested in the derivative of the band-limited signal represented by the samples.


Digital Distortion effect algorithm I read the DAFX book by Udo Zölzer about the distortion effect at pages 124-125
and it says that suitable simulation of the distortion is given by the function:

$$f(x)=\frac{x}{|x|}\left(1-e^{x^2/|x|}\right)$$

Can someone explain this formula and what kind of signal we get?

From what I understand 'x' is the sampled signal, so this is a sequence of numbers. What does |x| mean? Does it refer to the absolute value of x for each sampled value?

So if I want to implement this simulation of the distortion effect,

I need to know the length of x (It's given by the number of samples)
In a loop, I need to calculate this formula for each sample value
after the loop ends, I get the distorted signal (in a digital form)

After that, I need to convert it to an analog signal so I can hear it.

Thanks to the plot in Olli Niemitalo's answer I got convinced that the formula given in the book has a sign error. The non-linearity used for fuzz or distortion is always some type of smoothed clipping function, which compresses the input signal. So small input amplitudes experience little change whereas high input amplitudes are (more or less) softly clipped. And the figure shown in Olli's answer does exactly the opposite.

So I'm convinced that the correct formula should be

$$f(x)=\frac{x}{|x|}\left(1-e^{-x^2/|x|}\right)=
\text{sgn}(x)\left(1-e^{-|x|}\right)\tag{1}$$

For small values of $x$ we have $f(x)\approx\text{sgn}(x)|x|=x$, and for large (magnitude) values we get $f(x)\approx\text{sgn}(x)$, i.e., clipping.

This is a plot of the corrected non-linearity $f(x)$ (WolframAlpha):

The formula should also be simplified like the right-most expression in $(1)$, because a beginner might be inclined to literally implement the other formula and try to evaluate the terms $x/|x|$ and $x^2/|x|$, which is unnecessarily complex and also gives trouble when $x$ is close to zero. A typical implementation would look like this:

if (x > 0)
   y = 1 - exp(-x);
else
   y = -1 + exp(x);
end


Image cleaning - removing all regular shapes ( circles, squares, text, icons etc. ) I'd like an automatic method to clean images. I've been looking at the GIMP, and its command line (cli) interface, in scheme / lisp. It looks very promising, but I'm a complete novice, so I'd be grateful for advice.

I'm, essentially, wanting to remove all man-made artefacts from a large number of photos - which is why I need to automate it, it's easy enough with just the one.

I'd like to remove anything with a hard-edge (square, rectangle, or part of one), or with a smooth curve (circle, ellipse, or part of one), but not remove naturally curved things, like leaves, or fronds. Actually, since there may be bits of vaguely straight edge, or approximate arcs in the 'natural' figures, I can go further, and say that no straight edge, or arc, will be longer than .5 cm in the original - I hope that helps!

What's special about the areas removed is one of: There is a straight edged border (at least one). There is a border that is the arc of a circle - an interior or exterior border. The items left have irregular borders, with no straight edges and no arcs.

I've attached a before and after photo. In a perfect world, I'd like a script that takes the image on the left, and leaves the image on the right.

Any suggestions of the best way to go?

It'd be nice if the method could be one that can use a GPU.

I think the best way to get rid of those artifact shapes, is to implement an Shape Detection. Then, when you had coordinates of the unwanted shape, you can cut it out of your image. Then must fill its place with interpolating (through using background intensity values at the edges of the removed patch).

Now, there are various methods for shape detection and even OpenCV OpenCV Implementation  and Aforge  Aforge Implementation  support some of them. First you must get image edges, then using Hough transforms for lines and circles, detect shapes.

http://www.emgu.com/wiki/index.php/Shape_(Triangle,_Rectangle,_Circle,_Line)_Detection_in_CSharp

Finally as I said, cut out the image and fill out the removed patch through interpolation.


Is it possible for Matlab to read and decode an mp3 file? I have run the following Matlab code and I got a signal as shown in the Figure below.

[y, fs] = audioread(mp3file);
y1 = y(:,1);
y2 = y(:,2);
dt = 1/fs;
t1 = 0:dt:(length(y1)*dt)-dt;
t2 = 0:dt:(length(y2)*dt)-dt;
plot(t1,y1);
hold on
plot(t2,y2);

Is this a signal generate from a raw mp3 file or it is generated from a Matlab-decoded-mp3 file (which means Matlab has internally decoded the mp3 file)?

Is this a signal generate from a raw mp3 file or it is generated from a Matlab-decoded-mp3 file

The latter. $audioread()$ reads the raw MP3 bit stream, decodes it, and returns the decoded audio data as PCM data. If you want the raw bit stream use $fread()$ instead.


It is possible to obtain non-zero padding fft from zero padding fft? This is just a question out of my curiosity, normally we do not do this
I have a x(n)=[1 2 3]

y(n) = fft[1 2 3 0] =  6 + 0i  -2 - 2i   2 + 0i  -2 + 2i %zero padding
y'(n) = fft([1 2 3]) = 6 + 0i  -1.5 + 0.87i  -1.5 - 0.87i

I used Matlab to generate these results.
My question is there way to obtain y'(n) from y(n)?

Yes. Zero padding allows interpolation in the frequency domain. You can get the same result directly in the frequency domain using circular sin(x)/x interpolation. See for example
http://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula


How to analyze audio input/spectrum correctly I am quite new to Java and Signal Processing, but I was assigned a project dealing in audio processing. My topic given to me is a game that has a component of analyzing a song (any given song) and outputs data depending on the frequency currently being played (real time). ie: as the song plays, it outputs the current frequency (every second it outputs the current frequency).

I have been told on stackoverflow to use FFT. they say "just do a FFT" but that means nothing to me? HOW do you do an FFT? I have read tutorials, and basically understand what it is, but have no clue on how to implement it as in:

what is the data type of the audio file most suitable for such processing?
what exactly is the input of the FFT
how do I interpret the results

Can anybody suggest a simple to follow walk through/tutorial on how to process an audio signal? Additionally, if anybody knows a good implementation of FFT for Java, I would appreciate the suggestions.

Let's see. I have no idea about Java classes or what type of support it has for signal processing, but i will give you some guidelines. The Particular details of the implementation in the language, you have to figure out.

To do any kind of processing on an audio file you need the "raw" data, meaning an audio file wich has uncompressed audio samples. For exmaple, the WAV format, which normally has samples in 16bit signed 2's complement.

If you are going to use a computer and do offline processing, might be a good idea to cast that value to a double (Which in most processors today means 64bit floating point number).

Therefore you will have a continuos stream of double numbers coming from the file. Then you need to define a way to handle the continous data. A standard, very widely used way is to use a circular buffer (even if it is offline, i assume you want your code to be efficient, and loading the entire file to memory as an array is not the best solution). Either that or just a normal buffer, is your choice. The buffer length should be a power of 2 (Do to efficiency of the Cooley-Tukey radix-2 algorithm).

Now you need to do the actual FFT. This is merely a multiplicacion of the buffer (Which mathematically is a vector) by the FFT matrix. How this operation is actually performed in Java, no idea. In C it would be just to pass the pointer to the array and the length to an FFT routine that, either returns a pointer to dinamically allocated memory, or leaves the result in an array that you pass to it.

finally, you arrive at an array of M complex numbers (Assuming the length of the array/buffer/vector with the signal is M). And then you do whatever you want with it.

For instance, you could take the magnitude of each of the complex numbers, and find the maximum, to detect where the fundamental frequency could be (very approximate though).

Extras:
More advanced techniques to process would include a pre windowing to avoid leakeage, doing zero-padding to obtain more resolution on the windowed spectrum, etc.

Hope it Helps.


Help Understanding Convolution Filter Formulation I'm reading Computer Vision: Algorithms and Applications wich is available online as a PDF and on Chapter 3.2 page 111 it introduces convolution operators for images giving the formulation:

$g(i, j) = \Sigma_{k,l} f(i − k, j − l)h(k, l) = \Sigma_{k,l}f(k, l)h(i − k, j − l)$

since $(i,j)$ are pixel coordinates $f$ is the pixel intensity function for the image we're convolving and $h$ is the function for the convolution matrix it makes sense to me that $g(i,j)$ is equal to the first sum expression but I can't wrap my head around how this is equivalent to the second sum expression. If someone could explain this I would really appreciate it.

Te elaborate on nikie's comment that it is just a change of variables, consider the following.

Suppose that the original image is a single pixel of intensity $1$ at point $(0,0)$
and all other pixels have intensity $0$. Then,
the result of the convolution of this image is $h(i,j)$, that is, this
single pixel in the original image produces pixel intensity $h(i,j)$ at
the point $(i,j)$.  The actual image is a whole bunch of pixels, though, at
different points, and the superposition principle says that the resulting
pixel intensity $g(i,j)$ at point $(i,j)$ is the sum of of the pixel intensities
produced by the various pixels in the original image.
There are two ways of calculating $g(i,j)$.
Let's go through the pixels in the original image in order of their
distance from the point $(i,j)$.  Well, the point $(i-k,j-l)$ is at distance
$(k,l)$ from $(i,j)$ and the pixel intensity $f(i-k,j-l)$ there
produces intensity $f(i-k,j-l)h(k,l)$ at $(i,j)$.  This gives the first
formula
$$g(i,j) = \sum_k\sum_l f(i-k,j-l)h(k,l).$$
Alternatively, we can go through the pixels in natural order and ask
what intensity the pixel $f(k,l)$ produces at the point $(i,j)$. Well,
the distance is $(i-k,j-l)$ and so the intensity is $f(k,l)h(i-k,j-l)$
giving
$$g(i,j) = \sum_k\sum_l f(k,l)h(i-k,j-l).$$


Reconstructing time domain signal with Hanning window I am currently working on a frequency domain real-time application on a digital signal processor. Currently for one time frame of my algorithm I read in time domain data into a buffer, perform a Fourier transform with a Hanning window (win), which I generate using MATLAB with the following code:

win = 2*hanning(N,'periodic')/N;

where N is my transform length, (usually 2048 in my case). The window coefficients are then stored in a header file as single floating point data to be included into a C project.

In the frequency domain I execute my algorithm, then convert it back to the time domain to be outputted in a third buffer. I am certain the forward and reverse transforms work correctly having thoroughly tested them (and having used them for several months).

However when I output the data I use the following 'window' to reverse the effects (i.e. the energy level of the time domain signal to remove the 'oscillating' effect of the reconstructed signal) of the original window that I applied (in MATLAB code):

win_out = win ./ (win.^2);

I've seen the divide by energy (win.^2) in other peoples code, and assume that the line of MATLAB code above is equivalent to (it gives the same result of course):

win_out = 1./ win;

When I implement this in my real time application the values on the edge of my 'output' window (win_out) are very large (of the factor x10^8, see image below) and are in fact too big. This results in the edges of my time period being amplified too much resulting in a 'clippy' effect in my time domain signal when I listen to it.

I've illustrated the problem below

Above is the original reconstructed time domain signal, if this was to be outputted it would have an oscillating effect on the output, which I try and fix by applying an 'inverse' window.

This is the inverse window that generated by the above MATLAB code.

Time domain signal with my 'inverse' window (note: spikes at each end).

First spike zoomed in.

NB: I don't use overlapping (yet!)

My question is: how do a choose an 'inverse' (if that is the correct term?) window to ensure those spikes aren't present in fully reconstructed time domain signal with window?

Many thanks !

You should be using overlap.  That ought to take care of the "oscillating" and make it unnecessary to undo the Hanning window.

The effect you are seeing comes about because the window values at the edges get very close to (and may even be) zero.  The window edges won't line up perfectly with the signal after you've done your processing in the frequency domain, so dividing by those very small values will cause the peaks you see.

In your place, I would first make sure that I could do the windowed FFT/IFFT and recover the original signal.  When you can do that without getting the "oscillation" you can put your processing back in.


How is a PCM byte converted to a voltage level? 16-bit signed PCM would have possible values of ~-32,000 to ~+32,000. How is this value converted to a voltage level for the speaker?

If a DAC outputs 5VAC, the values would be ~+-5v. How is the binary value, such as 1101101011011101, converted to a voltage level?

There are quite a few different DAC topologies with different properties. The typical tradeoff is between resolution and speed and power usage -- higher speed means lower resolution, and higher resolution means higher power usage.

The simplest DAC is arguably the string DAC or Kelvin divider. If you consider the 3-bit DAC above, you can see that every possible value representable by 3 bits has a corresponding voltage level in the string. As indicated, you decode the 3-bit word to a switch a voltage level to the output, turning off all the others. For example, if you have the word 010, you turn on the 3rd switch, and if you have 101, you turn on the 6th.

You can extend this string to have $2^{16} = 65 536$ levels, but this make a very large circuit, and is not feasible in practice. The more compact R-2R ladder is then a better option, although a pure 16-bit R-2R ladders will have very high requirements for component tolerances. A combined string and R-2R topology is quite common for 16-bit.

For audio, the so-called $\Delta-\Sigma$ modulator is quite popular. In this case the output is switched, and averaging in the form of a low-pass filter is used to generate the desired levels.


Matlab - Constructing a Cosine wave with IFFT My apology if following question is too simplistic to be asked but I have the following code which generates a nice sine wave for me:

N=64;
X=zeros(N,1);
X(2)=-32j;
X(N)=32j;
x=ifft(X);
stem(real(x));

Now I was thinking I should be able to construct a cosine wave in the similar way by changing the X(2) and X(N) values into some real mirrored values say X(2)=-32 and X(N)=32. However, the output is not what I want i.e. a cosine wave

You must keep in mind that for a real-valued signal, second half of your spectrum is a complex conjugate of all values below the Nyquist frequency. In your case:

X(2) = -32
X(N) = 32

as you can see the coefficients are not the complex conjugate of each other. Because of that you are ending up with round-off errors, since two frequency components are cancelling each other out.

Here is some modification to your code:

clc, clear all, close all
fs = 32;
dt = 1/fs;

N=64;
X=zeros(N,1);
X(2)=32;
X(N)=conj(X(2));
s=ifft(X)';

td = linspace(0, (length(s)-1)*dt,     length(s));

stem(td, s, 'b', 'linewidth', 2)
grid on
legend({'reconstructed signal'})

Which gives:


MATLAB's $\tt cpsd$ and $\tt pwelch$ - different results for cross spectral power density I am wondering why the following code does not yield the same results for Sxy1 and Sxy2, where Sxy1 is calculated using the MATLAB's pwelch function, and Sxy2 is calculated using cpsd function.

fs = 1000;
t = 0:1/fs:5-1/fs;
win=hann(500);
overlap=300;
nfft=500;
fs=1000;

x = cos(2*pi*100*t)+randn(size(t));
y=sin(2*pi*160*t)+cos(2*pi*150*t)+randn(size(t));

[Sxx,f] = pwelch(x,win,overlap,nfft,fs);
[Syy,f] = pwelch(y,win,overlap,nfft,fs);

Sxy1=sqrt(Sxx.*Syy); %******

[Sxy2,f] = cpsd(x,y,win,overlap,nfft,fs);

Sxy2=abs(Sxy2); %******

Whilst I have not worked out how to calculate CPSD from power spectra, I have confirmed the following two methods are equivalent (using amplitude spectra, rather than power spectra or cross-correlation):

Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));
y=sin(2*pi*160*t)+cos(2*pi*150*t)+randn(size(t));
N = length(x);
nfft=500;
lfft=2*nfft;
overlap=0;
win=rectwin(N);

xdft1 = fft(x,2*nfft);  <outputs two-sided complex amplitude spectra
ydft1 = fft(y,2*nfft);
Sxy1a=xdft1.*conj(ydft1);
Sxy1b=Sxy1a./N^2;
Sxy1c=Sxy1b(1:lfft/2+1);
Sxy1c(2:end-1)=2*Sxy1c(2:end-1);
Sxy1=abs(Sxy1c)';

%CPSD Method 2
Sxy2a=cpsd(x,y,win,overlap,2*nfft,Fs);   < uses the cpsd function with Welch’s overlapped method
Sxy2=abs(Sxy2a);


Evaluating clustering result I am working with phoneme clustering. I have been using $k$-means clustering have done $k$-means on 2 different data: spectrogram vector and MFCCs. I would like to know how I can evaluate the result of clustering quality? I mean I need to know standard ways to validate the clustering. It would great if someone explains their experience.

One standard way to evaluate the quality of such techniques is to look at the distribution of correct clusterings versus incorrect clusterings.

This can be quantified by looking at the precision, recall and quality of the clusterings... but assumes that you have a "ground truth" (i.e. that you know the real cluster the data point belongs to).

Precision: $p = \frac{TP}{TP + FP}$, the percentage of positives that are true positives (and not mis-classified negatives).

Recall: $r = \frac{TP}{TP + FN}$, the percentage of true positives that were correctly identified.

Quality (AKA accuracy) = $q = \frac{TP + TN}{TP + TN + FP + FN}$, the percentage of correctly classified items from all classified items.

Here $TP$ is the number of true positives, $TN$ the number of true negatives, $FP$ the number of false positives, $FN$ is the number of false negatives.

If you have many clusters and need to evaluate their interaction, then a confusion matrix may be more useful.

If you want to go into this in even more detail, I'd suggest reading the paper:

  Andrew P. Bradley, The use of the area under the ROC curve in the evaluation of machine learning algorithms, Pattern Recognition, Volume 30, Issue 7, July 1997, Pages 1145-1159, ISSN 0031-3203, http://dx.doi.org/10.1016/S0031-3203(96)00142-2.
  (http://www.sciencedirect.com/science/article/pii/S0031320396001422)
  Keywords: The ROC curve; The area under the ROC curve (AUC); Accuracy measures; Cross-validation; Wilcoxon statistic; Standard error

A pre-publication copy of which is available.


properties of System I am studying Signals and Systems. Anyone, please help me to answer the problems.

$y(t)=dx(t)/dt$ is time invariant:  True or False ?
$h(t)$ is the impulse response of an LTI system. If $h(t)$ is periodic and non-zero, the system is unstable:  True or False?

For question 1: apply the definition of time invariant: find the output as normal; find the output with the same input but delayed by $T$
$$
y_1(t) = \frac{dx(t)}{dt}\\
y_2(t) = \frac{dx(t-T)}{dt}\\
$$
Does $y_2(t) = y_1(t-T)$ ?

For question 2: Find a definition of stability and apply it. For example, for a system to be BIBO stable it needs to have

$$
\int_{-\infty}^{+\infty} \left|h(t)\right| dt < \infty
$$

If $h(t)$ is periodic and non-trivial (zero), can that be true?


Why is my discretized transfer function unstable when my discretized state-space model is stable? I am working with a 30-state, 14-input linear model that is described by a state-space model:

model_state_space = ss(A, B, C, D);

The model is extremely slow (it describes the thermal behavior of a building), so when I discretize it I choose a sample time of 900 seconds and plot the step response from one of the inputs:

[y,t] = step(model_state_space);
plot(t, y(:,:,end))

I obtain a similar plot when I compute the step response of the discretized state-space model, or of the equivalent transfer function continuous model:

hold on
[y,t] = step(c2d(model_state_space, 900));
plot(t,y(:,:,end), 'r')
[y,t] = step(tf(model_state_space));
plot(t,y(:,:,end), 'k')

But when I compute the impulse response of the discretized transfer function, the system appears unstable:

[y,t] = step(c2d(tf(model_state_space), 900));
plot(t,y(:,:,end), 'g')

Can someone please explain why this is so?

You're using a high order model so you might be running into numerical precision issues.

Make sure the realization of model_state_space is minimal. You can do this with

min_model_state_space = minreal(model_state_space)

Also, not all state-space realizations are created equally. Some have better numerical properties, such as a balanced realization:

[balanced_model_state_space, G] = balreal(model_state_space)

If the vector G contains small values it means you can reduce the order of your model without giving up accuracy. Look up the matlab function balred.

Just from looking at your step response plot you probably don't need a 30th order model (1st or 2nd order might even be enough).


Reconstruction of voice from distorted recordings I would like to know how much of the sound of my voice getting through various types of walls can be reconstructed with easily obtainable software and algorithms.

My plan is to place a microphone on the other side of the walls I'm interested in and record the sound getting through. I would then apply reconstruction and try to understand my own speech. Since I am the one talking and I know what I was saying, this would be best case.

Are there ready algorithms, as standalone software, libraries or easy to implement in common software packages like Python, Matlab, Octave, ... (I'm talking about script-kiddie level of knowledge) that I can try to reconstruct offline the speech from the distorted recording?

This is a difficult problem. The wall itself is simply a linear time invariant filter. Typically it's just a low pass filter and can theoretically be compensated with a matching inverse filter.

However in practice the attenuation at most frequencies is quite severe. Unless the recording room is extremely quiet and the recording extremely low noise, the remaining speech signal will be below the noise floor and hence impossible to recover.


Finding a good inverse for an ill-conditioned matrix transformation I have a time-series observation dataset that has been distorted. I want to recover the best approximation of the original signal as possible. Disclaimer:: I know only the basics of linear algebra, so please bear with me.

I have a good model of the distortion, represented by a square matrix. In theory, if I can find an inverse of the transformation, I can recover the original signal. However, the model of the distortion as a matrix is ill-conditioned (i.e. almost singular).

Is it possible to take this matrix model, and generate an invertible matrix that is an approximation of the original matrix?

The usual way to approach inverting a matrix that is rank deficient is to use a generalized inverse.

This means if your matrix equation is:
$$
y = \mathbf{A}x
$$
where $\mathbf{A}$ is an $n \times m$ matrix, then you can use, for example, the Moore-Penrose pseudo-inverse, $\mathbf{A}^{+}$:
$$
\mathbf{A}^{+} = (\mathbf{A}^*\mathbf{A})^{-1} \mathbf{A}^*
$$
if a left inverse is required or
$$
\mathbf{A}^{+} = \mathbf{A}^*(\mathbf{A}\mathbf{A}^*)^{-1}
$$
for the right inverse.

If $n=m$ then either be used.

However, as per the link above, there are many different possibilities for selecting an inverse in the rank-deficient / non-square case because the system of equations is underdetermined.

As per @MBaz's comment, you can calculate this using the singular value decomposition of $\mathbf{A}$.


how to correct minBLEP sawtooth oscillator DC bias In my implementation of minBLEP, overlapping minBLEPs result in a DC bias as the frequency increases. This is because each minBLEP correction (i.e. the minBLEP minus the aliased step) has a DC bias.

Here's a frequency sweep from 100Hz to 15kHz. The second signal shows the number of active BLEPs.

Interestingly, the minBLEP paper (http://www.cs.cmu.edu/~eli/papers/icmc01-hardsync.pdf) acts as if the DC woes are all solved by switching from integrating impulses to placing BLEPs. There's no mention of correcting DC when lots of BLEPs overlap.

Should I just add a DC-blocker at the end? Is that what is generally done?

Also, why is it important that the BLEPs are minimum phase? Using linear phase BLEPs would alleviate the DC bias, at the cost of delaying the oscillator output by half the length of the BLEP.

If your oscillator is steadily accumulating DC over time, adding a DC blocker at the output may not be an elegant solution since your oscillator will keep getting more and more DC.  For floating point this may or may not be a huge problem (I wouldn't do it), but for fixed point you will eventually start clipping or wrapping.

If the oscillator has a frequency dependent DC offset that doesn't creep up over time, a DC blocker is probably a good solution.  If you've implemented the BLEPs correctly, this should be the case.

The original reason given for BLEPs to be minimum phase was to implement oscillator sync.  To implement this correctly with linear phase BLEPs, you should use lookahead so you can start the BLEP before the aliased reset occurs.

However, apparently this author would rather not have any lookahead for some reason and therefore designed minBLEPs so that you would have to lookahead less far into the future.  If I remember correctly, in the paper they say that it's okay to just forget about lookahead if you use minimum phase, but obviously this isn't perfectly true as minimum phase filters still have some delay.  You should experiment yourself to see if the approximation works for you.

In my opinion, it's probably a good idea to use linear phase BLEPs and lookahead appropriately for sync.  We're talking a few samples of latency here for reasonable BLEPs.


What's the right scale to represent a frequency spectrum? To represent a 16 bits temporal spectrum, it's "easy" : my scale is going from -32000 to +32000. But What's the scale for a frequency spectrum ?

I made that try : I calculated a pure single frequency sine signal and calculated it's energy for a period. So, I supposed that any signal generated on such a period - without saturation, of course - could have more energy. So, I considered the ratio

  20.log10 (signal_energy/maximum_theorical_energy)

as my 0 dbFS.

Then, I made few experimentations, generating signals for different kind of frequencies and, when I had a single sine, each time, I reached the 0dbFS.

But, of course, a music have thousands of frequencies. So, when I represent a spectrum of a "music sound" signal, I rarely reach more than -30 db and my signal is more or less around -60 db for each frequency. It seems very low.

So, my question is : is it correct to represent a signal with such a measurement ? If not, how do represent the "values" for each frequencies on the graph ? With which scale ?

Thank you for all your answers.

FFT scaling for audio signals is indeed difficult. Sometimes the energy is distributed over many frequencies, sometimes the energy is only at few or just one frequency. This can result in a very high "crest factor". i.e the ratio of the peak value to the average power (RMS) value.

This mainly depends on the the length of the FFT: the longer the FFT, the larger the crest factor will get.

The best ways to deal with this is:

Use floating point representation and math if any possible
Choose an FFT length that is optimal for your application and not longer than you really need

Other than that, the differences that you see are real. FFT is a metric of "spectral density", and for pure sine waves this approaches infinity, while it is quite low for noise-like signals.


What is the step response curve of a second order low pass filter? And what if it is set for resonance? I am doing a lot of audio synthesis DSP where I need to use low pass filters to shape the decay of an impulse. I understand that a one pole filter fed a discontinuity, say from a steady state of 1's to a steady state of 0's, will decay with a perfectly exponential curve. I believe the decay will be such that the signal will reach an amplitude of $1/e$ at $time = 1/(2*pi*f)$.

This makes conceptualizing one pole low pass impulse modeling very easy. But I am trying to understand how other filter orders will do the same.

I understand that a second order filter can be approximated running two one pole filters in sequence (ie. each will attenuate by 6 dB/oct, leading to a 12 dB/oct curve). If that is the case, then I would expect the step response of a second order filter without resonance to behave the same as a one pole filter and remain exponential, only that much faster.

So perhaps the time to 1/e would be:

$time = 1/(2*pi*f)^2$

Would that be correct?

I asked here about a lower slope filter and no one has replied so I am guessing it is not common knowledge:
A one-pole LPF (6 dB/oct) has a step response to $1/e$ amplitude of $time = 1/(2∗pi∗f)$. What would the response time be of a 3 dB/oct filter?

But I presume the same principle would then apply. A 3 dB/oct low pass filter would have a time constant of:

$time = 1/(2*pi*f)^{0.5}$

And if this is all correct, then we could say the time constants of any non-resonant low pass filter would be roughly:

$time = 1/(2*pi*f)^{filter-order}$

What do you think? Does all this sound correct? If so, there is no difference between low pass filtering with any order of non-resonant filter, as they can all be set up to create the same outcome with different parameters.

Lastly, then the remaining question would be: What does the step response of a resonant second order low pass filter look like, again going from steady 1's to steady 0's? I presume the resonance ruins the precisely exponential decay, but in what manner? Does it create a more compressed curve? Will it resonate creating a wobbling of the output? Will it dip below zero as it tries to settle out?

I tried testing it with a resonant second order LPF but I was just getting almost random maximum amplitudes coming out. Very unpredictable. I'm not sure if I did something wrong or that's to be expected from the resonance.

Thanks for any help understanding all this.

I can help with some of your multiple questions. First, a cascade of n buffered RC low pass filters (LPFs), a so-called n-th order synchronous LPF, has the impulse response and step response shown in the screenshot below:

This is a screenshot from my paper 1 referenced at the bottom. All R values are the same and all C values are the same. Buffering simply means there is a unity gain buffer between each filter and on the first input and final output as well. So no filter is loaded or gets loaded. Notice that the conventional unit step function goes from 0 to 1, so if you want the response to be 1 to 0, just subtract y(t) from 1. Then you can find 1/e times as a function of n, etc.

With regard to the four "time = ..." equations you give, only the first one, for n = 1, is correct. In general, second (and higher) order filters do not all behave the same. The simple RC LPF is the first (primitive) member of several filter families. But things get much more complicated (and useful) as the various properties of the filter families are taken into account and utilized.

Note that a cascade of RC LPFs cannot have overshoot or undershoot. On the other hand, Butterworth LPFs can have "ring" in response to a unit step input. This is nicely illustrated by the following figure from Blinchikoff and Zverev 2:

References:

1  E. Voigtman, J.D. Winefordner, “Low-pass filters for signal averaging”, Rev. Sci. Instrum. 57 (1986) 957-966.

2  H.J. Blinchikoff, A.I. Zverev, "Filtering in the Time and Frequency Domains", Wiley-Interscience, John Wiley & Sons, NY, ©1976, p. 114.


Difficulty regarding grasping concept of fft butterfly? The  term "butterfly " is often used in association with fft texts
But what exactly is a butterfly
?a matrix/group of i/o's is known as butterfly? Such as 2 input 2 output butterflies in stage 1 of attached photo
How we can differentiate between butterflies of stage 1 and stage 2?other than number of i/o's?

A FFT "butterfly" is the name for an algorithmic structure inside the FFT. It

has two complex inputs
two complex outputs
one complex multiply
sum and difference of two complex numbers

There are two basic types. The decimation-in-time butterfly does the multiply first
$$y_0 = x_0 + x_1 \cdot W$$
$$ y_1 = x_0 - x_1 \cdot W $$
Decimation-in-Frequency computes the sum/difference operation first
$$ y_0 = (x_0 + x_1) $$
$$ y_1 = (x_0 - x_1) \cdot W $$
It's useful to write this in matrix notation. Decimation-in-time:
$$ \begin{pmatrix}
y_0 \\ y_1
\end{pmatrix} = \begin{pmatrix}
1 & W \\ 1 & -W
\end{pmatrix} \cdot
\begin{pmatrix}
x_0 \\ x_1
\end{pmatrix} $$
and decimation-in-frequency:
$$\begin{pmatrix}
y_0 \\ y_1
\end{pmatrix} = \begin{pmatrix}
1 & 1 \\ W & -W
\end{pmatrix} \cdot
\begin{pmatrix}
x_0 \\ x_1
\end{pmatrix} $$
So the matrices are transposed of each other.
Update:
The twiddle factor $W$ is a function of which stage you are in and which butterfly inside the stage. For decimation-in-time, $W$ simply alternates between $+1$ and $-1$, i.e. $W = [+1,-1,+1,-1 ...]$ so there is no need for an actual multiplication. Similarly for the second stage we have $W=[1, -j, -1, j, 1, -j, ...]$ which also doesn't require multiplication. This property can be used to further optimize the implementation and so in many Implementation these stages are hand-coded


Determine constant $A$ such that $x[n] = x[n] \star x[n]$ Let $x[n] = A\delta[n] - \frac{\sin(\frac{3n}{2})}{\pi n}$. Determine constant $A$ such that for all $n$ $$x[n] = x[n] \star x[n] \tag{1}$$
I think it's not possible since $(1)$ leads to $$X(e^{j\omega}) = X(e^{j\omega})X(e^{j\omega})$$ And this means $X(e^{j\omega}) = 1$ or $X(e^{j\omega}) = 0$. Also $$X(e^{j\omega}) =  \begin{cases} A - 1   &0\le | \omega| \le \frac{3}{2}   \\ A &   \frac{3}{2}\lt | \omega| \le \pi \end{cases}$$
It means no value of $A$ works. I don't know whether is my answer correct. Maybe I've neglected something.

Systematically, you could just solve the problem by finding a solution $A$ that satisfies the following two equations:
$$(A-1)^2=A-1\\A^2=A$$


Why do we use "error quaternions" for control of attitude? I'm working on an attitude control system in terms of quaternions, and i'm basing my approach on the book Space Vehicle Dynamics and Control (2nd ed.).  If we consider a proportionate controller, i.e.
$$\mathbf{u} = -\mathbf{K}_p \mathbf{q}_e $$
it would seem logical to me to define the error quaternions to be the first three elements of
$$\begin{bmatrix}
  q_{1e} \\
  q_{2e} \\
  q_{3e} \\
  q_{4e} \\
\end{bmatrix} = \begin{bmatrix}
  q_{1} \\
  q_{2} \\
  q_{3} \\
  q_{4} \\
\end{bmatrix} - \begin{bmatrix}
  q_{1c} \\
  q_{2c} \\
  q_{3c} \\
  q_{4c} \\
\end{bmatrix}$$
Where $q_i$ are the actual quaternions and $q_c$ are the desired (commanded) quaternions. Yet in the book, they instead define
$$\begin{bmatrix}
  q_{1e} \\
  q_{2e} \\
  q_{3e} \\
  q_{4e} \\
\end{bmatrix} = \begin{bmatrix}
  q_{4c} & q_{3c} & -q_{2c} & -q_{1c}  \\
  -q_{3c} & q_{4c} & q_{1c} & -q_{2c} \\
  q_{2c} & -q_{1c} & q_{4c} & -q_{3c} \\
  q_{1c} & q_{2c} & q_{3c} & q_{4c}\\
\end{bmatrix}\begin{bmatrix}
  q_{1} \\
  q_{2} \\
  q_{3} \\
  q_{4} \\
\end{bmatrix}$$
To me it seems like we do this because there are multiple possible quaternions that correspond to the same orientation, but I have not yet found a clear explanation for this. Can anyone enlighten me?

A possible explanation is that a set of 4 quaternions is not independent, so by converting them to error quaternions, we get 3 independent values that we can use to calculate 3 independent control actions (one for each axis).


What is the frequency transfer function when $x(t) = e^{-3t}u(t)$ and $y(t) = 2u(t)[e^{-t}-e^{-4t}]$ What I did was take the fourier transform of both $x(t)$ and $y(t)$ and then divided $Y(j\omega)/X(j\omega)$. So
$$Y(j\omega) = \frac{2}{1+j\omega}-\frac{2}{4+j\omega}\quad\text{and}\quad X(j\omega) = \frac{1}{3+j\omega}$$

However, the answer is of the form: $\displaystyle \frac{A_1}{1+j\omega}+\frac{A_2}{4+j\omega}$ where the $A$'s are constant. I am not getting the solution of this form. I got $Y(j\omega)$ in this form, but not $H(j\omega)$.

Am I doing something wrong?

Your answer is correct; it's just in a different form than the given answer. In order to get your answer in the given form, do the following:

rewrite $Y(j\omega)$ by combining the two terms: $Y(j\omega)=\displaystyle\frac{N(j\omega)}{(j\omega+1)(j\omega+4)}$, where $N(j\omega)$ is a (very simple) polynomial. I'm sure you know how to obtain $N(j\omega)$.
write the frequency response as $\displaystyle H(j\omega)=\frac{N(j\omega)(j\omega+3)}{(j\omega+1)(j\omega+4)}$ and use partial fraction expansion to obtain the answer in the given form.

Alternatively, you could also determine $h(t)$ via inverse Fourier transform of $H(j\omega)$ as you obtained it (and noting that multiplication by $j\omega$ corresponds to differentiation in the time domain). You will see that $h(t)$ is in the form $A_1e^{-t}u(t)+A_2e^{-4t}u(t)$, which directly leads to $H(j\omega)$ in the given form.


Why noise canceling headphones do not attenuate the signal itself? Noise canceling headphones use a microphone to capture noise in order to cancel it. How is it guaranteed that the microphone will not record the actual signal and attenuate it?

The microphone will record the noise on the outside and the speaker sits on the inside with a bit of passive sound shielding in between. Outside noise will reach the microphone unattenuated, but the actual signal will be attenuated by the passive shielding. The noise that reaches the speakers attenuated by the passive sound shielding and the circuit of the active noise canceling will simulate the effect of the passive shielding to generate the right amount of anti-noise to add to the actual signal for playback on the speaker. The circuit can also simulate how much of the actual signal reaches the microphone and compensate for that or it can just ignore that issue because this effect would only cause a small attenuation of actual signal.

Example:
Lets say our passive shielding reduces sound levels of signals traveling through it by 75%, so 25% reaches the other side. Actual signal on the inside without noise canceling would be 1.0, noise on the outside would be 1.0 as well. Without noise canceling noise level on the inside would 0.25 and actual signal on the outside would 0.25. So the microphone would record 1.0 noise + 0.25 noise, now when we turn on active noise canceling, as a first order model, the circuit would simulate the attenuation of the shielding and subtract 0.25 noise and 0.0625 actual signal. If the attenuation is simulated perfectly, on the inside we would get 0.9375 actual signal and no noise. SNR is significantly improved. Of course this is just a first order model, the attenuation provided by the passive shielding does not have a flat frequency response and the microphone has also some electrical noise, that will be added to the speaker signal. If the simulation of the passive attenuate is not perfect, we will also add extra noise. We also need to adjust for the frequency response and gain of microphone and speaker. Small errors can easily cause us to add more noise to the signal or subtract to much from it.