“Subtracting” one audio stream from another audio stream I am a long time software engineer but have practically zero experience with signal/audio processing. I am interested in learning about signal processing via a use-case we have for one of our audio components. NOTE: this is just a learning exercise for me... it is not a priority that the end result be useful.

We have a component (A) which produces audio (via a speaker). We have another component (B) which records that audio (via a simple microphone).

What I would like to do is use A to record B's audio. Then I would like to, if even possible, compare the two streams. The goal would be to remove the audio that was present in stream A - leaving the ambient conditions which existed during the recording - I understand full removal is not possible.

I realize that there are phase and magnitude issues. I also realize that it isn't just a simple matter of "subtracting" B from A. That said, my assumption is you can subtract A from A.

I would like to understand how to approach the problem. Again, this is a learning experience for me (there are no deadlines); I am more than willing to start from the beginning.

Any advice/suggestions would be much appreciated.

This setup shares some similarities with system identification problems, where $A$ would be the input of the LTI system you want to estimate the transfer function of, $B$ being the output; and the "ambient sound" being the additive noise. The LTI assumption is reasonable provided your converters/amplifiers/transducers are of decent quality.

So the steps would be:

Use a system identification technique to find the FIR filter $\hat{h}$ that minimizes the mean-square error between $\hat{h} \star A$ and $B$. A simple method, which might not be the most suitable here, is to divide the cross-correlation of $A$ and $B$ by the autocorrelation of $A$. Explanation here. The limitation is that it will not work well for long recordings (you might better compute your estimates on shorter segments and average them) - and that music is not the best "probe" signal to send into a system to estimate its response.
You can now use $\hat{h} \star A$ as an estimate of the original signal A as "heard" by the microphone and subtract it from $B$ to retrieve the ambient sound.

I gave a shot at this using a music clip (A), applying a reverb and slight amp model to simulate a speaker in a room, then mixing in a cat audio sample to get (B), then estimating an impulse response from the (A, B) pair, then subtracting the filtered A from B. This shows some results but a better FIR estimation technique might help here! (note that I truncated the estimated IR to its first 5000 samples to speed up computations).

Note that there are algorithms for doing this adaptively (such as LMS). This might be more suitable for your problem if $A$ and $B$ are processed in realtime rather than offline. Such algorithms form the basis of echo cancellation systems used in telecommunications.


the demonstration that states that the FT of the ACF function is the square of the DTFT of the signal I am following the book The Intuitive Guide to Fourier Analysis & Spectral Estimation with MATLAB. I am trying to selflearn the fourier analysis in matlab.
I got lost in one passage in the demonstration that states that the FT of the ACF function is the square of the DTFT of the signal. I have attached it here As you can see in the passage that I named 1. a delta_tau is missed in my opinion. Can you confirm that? More important in the point 2. in my opinion there should be a tau and not t. So I don't know how it is demonstrated this formula, because if in the passage 2. there is t the demonstration can be ended. I hope you can help me

You are right, the derivation is full of typos. The first equation below Eq. $(8.39)$ should read

$$\int_{-\infty}^{\infty}x(t+\tau)e^{\color{red}{-}j\omega\tau}d\tau=X(\omega)e^{j\omega \color{red}{t}}\tag{1}$$

Substituting into $(8.39)$ gives

$$\begin{align}\mathcal{F}\big\{R(\tau)\big\}&=\int_{-\infty}^{\infty}x(t)X(\omega)e^{j\omega t}dt\\&=X(\omega)X(-\omega)=|X(\omega)|^2\end{align}\tag{2}$$

where the last equality is only true for real-valued $x(t)$. However, the overall result is also true for complex-valued $x(t)$ because in that case the ACF is defined differently:

$$\mathcal{F}\big\{R(\tau)\big\}=\int_{-\infty}^{\infty}x^*(t)x(t+\tau)dt\tag{3}$$


Necessary image processing operations before template matching Well, this may show I didn't take image processing class.

I have this template to match an object in images. But real object in images, for example, does not look like that but the second image. What operations do I need to take for the images like the second one to be ready for template matching? FYI, the third image shows a ideal matching.

In another idea, I fit the template to the image. That is, I give four corner coordinates on the image to match those in the template. Then it will compute a way to fit the template to the area.

UPDATE: Problem solved. Thanks to the answer I selected below. I successfully managed to apply my template on the image. The last image shows the result.

The sample image you posted has relatively strong perspective (it is not imaged straight from the direction of surface normal) which can cause problems with template matching techniques witch use block processing. I assume that you have to take the image with strong perspective so first thing we want to do is estimate image transformation which will remove the perspective projection and will result rectified image ("straight image").

Correcting perspective projection when shape is known

The ultimate goal of the preprocessing is to find corners of that gray rectangle which is the background the color chart. I started by simply thresholding with constant value. Literature/web has lots of information on how to do thresholding when the image has uneven lighting, but if you can, try to correct the imaging setup such that the lighting is relatively even.

I assumed that the color chart takes quite large portion of the screen, so the it is probably biggest continuous region of the image. There exist many black and white labeling functions (MATLABs bwlabel, IPPs LabelMarkers, OpenCV has external library cvblob) which will give each distinct connected region its own index. After the labeling you can calculate the areas of the connected components using a histogram and then choose the component with largest area.

It is usually good idea to fill the components such that component does not have holes, this reduces the number of edge pixels.

Now we can use edge finding operators (I used gradient based method, but you can also trace the edge by following white pixels that are connected to black pixels) to get pixels belonging to the boarder.

As the object we are interested is rectangle, I fitted four lines to the edge pixels using robust line fitting (RANSAC). After finding new line, I removed all the points that were close to the found line, this forces subsequent fitting operations to return other edges of the object. Finally when all the edges were found, I computed all possible intersections between the set of lines. From the set of resulting points I removed all points outside of the image.

Now we have the corners of the color chart and as we know the dimensions of the chart, we can calculate the transformation matrix $H$ (homography, projection between two planes. use direct linear transformation (DLT) to solve $H$) between the millimeter coordinates and the pixel coordinates.

$\textbf{x} = H\textbf{X}$

This transformation matrix can be then used to transform the original image such that color chart is imaged straight from the direction of the surface normal.

As can be seen only radial distortion is left. When we used transformation $H$ for image transformation (2D interpolation), we chose the interpolation coordinates such, that the we know the precise resolution of the image, meaning that we now know for example the size of the letters.

If we still want to do some template matching, we can use any decent template matching method. The method does not have to be rotation/scale invariant as the image is already rectified.


Correlating a sequence or a code with itself at QPSK symbol rate I am sending a QPSK signal that represents a certain sequence or code and  that needs to be correlated with itself on the receiver (there are other sequences/codes in the system, and they are all orthogonal). The result of the correlation should  indicate which code was sent . Suppose an AWGN channel, and at the output of the I and Q matched filters on the receiver  I have one sample/symbol that is a 'sufficient statistic" , so one sample on I and one sample on Q , and suppose timing synchronization is perfect and no other impairments besides AWGN (SNR is high enough), my question is the following:  Can correlation take place at the symbol rate ? or does one need more than one complex sample/ symbol in order to do the correlation?

The answer is yes, absolutely 1 sample per symbol is sufficient if your symbol timing and frequency offsets are sufficiently accurate (they need not be perfect).

Regarding frequency accuracy: The magnitude of correlation vs frequency offset is a Sinc function with first null at 1/T where T is the length of the sequence being correlated.

I started to develop this point in another recent post (Derivation of the Optimal Matched Filter - Convolution vs. Correlation) showing the first three samples of a correlation when a frequency offset exists. The results in this post showing correlation without a frequency offset should make it very clear why 1 sample (the correct sample) per symbol is sufficient.  In the plot we see the resulting Sinc response of correlation vs frequency offset with the first null at 1/T = 1 Hz. If we needed a wider tolerance to frequency offset, we could easily correlation over 100 samples in the case to then push the null out to 1/.1 = 10 Hz. Of course we would lose 10 dB in processing gain (see prior post referenced) so we see the trade of SNR needed versus the frequency inaccuracy that can be tolerated.

This next graphic continues that example to show the case of a correlation over 1000 samples vs frequency offset when the sampling rate is 1 KHz. This means that T = 1 second in this case, and we see that if we considered correlating a 3 Hz signal to a 4 Hz signal, that after summing 1000 products, the summation will be back at the origin (correlation null).

The exact relationship for the correlation versus timing offset depends very much on the pulse shape of each chip in your signal (the change in magnitude for each sample due to timing offset), and the concern therefore would be both with a timing offset in time, and a timing offset in frequency (this is different from the frequency offset above in that an offset in frequency of the timing clock will cause it's sample position within the symbol to roll, while a frequency offset in the carrier will cause the complex IQ correlation to rotate on the IQ plane.).

Both cases are easily computed by knowing that correlation is a sum of products, and correlation sequences typically have the same number of 1's and 0's (or close to the same number, often they are an odd number so in that case would have one greater 1 than 0 or vice versa).


Is it possible to use a Volterra series to generate subharmonics? Taylor series can create harmonics of the frequencies in an input signal. I'm wondering if it is likewise possible to use Volterra series to create sub-harmonics of the frequencies in that signal.

Some degree of finesse is required in posing the question:

Taylor series don't only create harmonics; there's intermodulation distortion as well. Likewise, it's ok with me if a suitable Volterra series doesn't only create sub-harmonics, but also creates something similar to intermodulation distortion, or even some degree of harmonic distortion.

However, I would like to exclude trivial cases where the input signal gets so wrecked that it ends up looking like white noise, which would obviously create sub-harmonics (along with every other frequency you can imagine!). I'm not sure how to formalize this requirement precisely, other than to say that I hope it's clear what I'm driving at.

Here's one way to formalize the behaviour I want:

Given an input $\cos(\omega t)$, yields an output containing $\cos(\frac{\omega}{n} t)$ for
some fixed $n$, ideally a natural number
Given a general input, a sub-harmonic is generated corresponding to every frequency in the input
Other "intermodulation distortion"-type artifacts are allowed, "within reason"

The answer can be expressed for either discrete or continuous signals.

Check out this paper.  I would have made a comment but not high enough rep.

http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=1211087&url=http%3A%2F%2Fieeexplore.ieee.org%2Fiel5%2F81%2F27258%2F01211087

Looks like you need multiple in to get subharmonics in Volterra series

The abstract states "Subharmonic generation is a complex nonlinear phenomenon which can arise from nonlinear oscillations, bifurcation and chaos. It is well known that single-input-single-output Volterra series cannot currently be applied to model systems which exhibit subharmonics. A new modeling alternative is introduced in this paper which overcomes these restrictions by using local multiple input single output Volterra models. The generalized frequency-response functions can then be applied to interpret systems with subharmonics in the frequency domain."


How to Detect Gradients and Edges in Images? I want to be able to find points in images that are the centre of a radial gradient like the one shown in the left picture below. Any ideas on how I could use a Hough transform or some other computer vision method?

Thanks

example search image:

I was working in opencv and trying to find the peak of a gradient generated by distance transform. I realised that using morphologic operations (erosion/dilatation) in grey-scsale images was very useful in this case.
If you erode dilate a grey-scale image, any pixel is going to take the value of the lower/highest neighbour. You can therefore find peaks of intensity in gradients by subtracting grey-scale image from the same dilated/eroded image. Here is my result:

And a way to do it in OpenCV/Cpp:

#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"

int main( int argc, char** argv ){

    cv::Mat objects, img ,peaks,BGR;
    std::vector<std::vector<cv::Point> > contours;
    /* Reads the image*/
    BGR=cv::imread(argv[1]);
    /* Converts it to Grayscale*/
    cv::cvtColor(BGR,img,CV_BGR2GRAY);
    /* Devine where are the objects*/
    cv::threshold(img,objects,0,255,cv::THRESH_BINARY);
    /* In order to find the local maxima, "distance"
     * is subtracted from the result of the dilatation of
     * "distance". All the peaks keep the save value */
    cv::dilate(img,peaks,cv::Mat(),cv::Point(-1,-1),3);
    cv::dilate(objects,objects,cv::Mat(),cv::Point(-1,-1),3);

    /* Now all the peaks should be exactely 0*/
    peaks=peaks-img;

    /* And the non-peaks 255*/
    cv::threshold(peaks,peaks,0,255,cv::THRESH_BINARY);
    peaks.convertTo(peaks,CV_8U);

    /* Only the zero values of "peaks" that are non-zero
     * in "objects" are the real peaks*/
    cv::bitwise_xor(peaks,objects,peaks);

    /* The peaks that are distant from less than
     * 2 pixels are merged by dilatation */
    cv::dilate(peaks,peaks,cv::Mat(),cv::Point(-1,-1),1);

    /* In order to map the peaks, findContours() is used.
     * The results are stored in "contours" */
    cv::findContours(peaks, contours, CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE);
    /* just draw them and save the image */
    cv::drawContours(BGR,contours,-1,cv::Scalar(255,0,0),-1);
    cv::imwrite("result.png",BGR);

    return 1;
}


Effect of redundant training data in HMM-based speech recognizer/synthesizer? How are redundant training data handled during the training stage?

For example, assume we have one observation for phone $\theta$ in the training set.

Then the training (for a monophone) is done with:

$$\lambda_{max}^\theta = \text{arg} \max_\lambda p(\bf{O}|\lambda) $$

where $\bf{O}$ is the observation sequence and $\lambda$ is the HMM model.

This is straightforward when only a single observation is available from the training set. But what happens when there are multiple observations for phone $\theta$ in the training set? How is $\bf{O}$ adjusted?

I guess you confuse GMM and HMM trainings. Although in both cases EM algorithm is employed, Baum-Welch is used for HMM training.


Timing error detector used for symbol timing synchronization The Gardner timing error detector is independent of carrier phase and frequency error, and is also quite immune to self noise (to my understanding, I may wrong about that). Can one assume that it is a recommended  timing error detector for these reasons?
Please if the answer is going to be "it depends on your system" with no further brief input/explanation , then abstain please.

It depends on your system :)

Seriously, it depends. Gardner is a pretty solid algorithm with a simple implementation, especially if you're using BPSK or QPSK. It can even be used in some MIMO systems. Its main disadvantage is that it requires two samples per symbol. M&M requires only one, but it has its own disadvantages.

I think that these days the preferred solution is some variation of the polyphase detectors first proposed by fred harris. However, if you're experimenting and building your own system and you want to keep it simple, Gardner is a safe bet.


Beginner - use FFT to display spectrum/equalizer At first, I want to say that I'm totally beginner in DSP and I don't have enough math education to understand Fast-Fourier transformation.

I'm a programmer and I want to archieve something like this:

My example file is in 8kHz, 16-bit samples, but currently, I'm trying to do with my own raw data - one sin(x) function. I know that I need to pass N samples to FFT function, but how interpret result? As I see, FFT need some input data (N samples) and gives N output data.

How I can know where is frequency and scale it - for example - to 5 scopes like 0 - 200Hz, 200 - 400Hz etc.

Or maybe I don't understand FFT correctly? I think that I can pass some samples and get information  "how much" is each frequency.

Technically, I'm using FFTW3 library.

I think we need a little more information to fully answer your question.  Most FFT routing (like FFTW) provide many different ways to call the transform, and the return result will depend on exactly what routine you call, and how you call it.

Having said the standard disclaimer, if you perform an FFT on $N$ real data points, you will get back $N/2$ unique complex points.  The basic algorithm though will assume your $N$ real points are $N$ complex points with zero imaginary parts.  Many routines will thus return $N$ complex points, but the second half will be (in the case of real input) hermitian symmetric (that is, point $N/2+i$ will be the complex conjugate of point $N/2-j$).

So, start by recognizing that your result is complex numbers, not real, and you will need to at least plot the "square" ($xx^*$) or each point.  If you have the real and imaginary points returned as $X_r$ and $X_i$, then you want to plot $X_r^2+X_i^2$.

So you are probably pretty close.  You will want to consider a few other things: windowing your input data, and averaging your output data.  To do the first, you need to perform the transform on weighted data, rather than your direct data.  See http://en.wikipedia.org/wiki/Window_function for more info.  You will simply hand the FFT function data of the form $W_k   {x_k}$ instead of just the $x_k$ (where $k$ is the index number).  This will be particularly important if you use a fractional frequency in your sin function (try it, you'll see).

The averaging step involves averaging the squared output. You will simply want to average a few outputs, point by point, together.  You will also often want to plot the output as a logarithm though, instead of simply plotting the squared output.


Interpolation of a gray scale image in matlab I want to fill the black areas with the value of neighboring pixels. Kind of interpolation. Can any one suggest me how Can I do in matlab.

Try the imfill function, see http://www.mathworks.com/help/toolbox/images/ref/imfill.html


Fourier Transform Computer Vision Textbook I am required extremely fast to fill the gap and learn Fourier Transform with application in Computer Vision.

It was easy to find mathematical aspects of Fourier Transform, but I am more interested in application in Computer Vision for a newbie with a lot of examples.

If you familiar with a textbook or have any other source of information about Fourier Transform well explained with examples in the field of Computer Vision, please share it with us.

Thats funny,  I was exactly trying to fill the same gap.
I found the several sites which helped me

this is one of them :
http://homepages.inf.ed.ac.uk/rbf/HIPR2/fourier.htm

Moreover , I recommend on creating a black image and playing with it.

e.g.

black = false(500,500);
black(10:40:end,:) = true;
black = imrotate(black , 45 , 'crop');
ffim = fftshift ( abs ( fft2 ( black ) ) );
lffim = log ( ffim ) ;
fs(lffim,[]);


Generate 3D image from Multiple 2D Images If we're to capture an image of a body using four (4) WebCam with a 16 megapixel (1 in the front, 1 in the back and 2 sides(left and right)). Is it possible to make a 3D out of it using image processing?

You can make a 3D model of the scene (<=> the depth) if the two cameras are looking to the same scene. After that you detect/extract the commun points between the two images and do some triangulation so you can recover the z component.
From your installation as you discribe it you can only enlarge the field of view


DFT dirac simplification I found this equality in the middle of my DFT calculation. Can anyone help me to proof this.

$ \frac{1}{N} \sum\limits_{k=0}^{N-1} \left(-W_N^n \right)^k \ = \ \delta[n-0.5 N] $

where

$ W_N = e^{j 2 \pi/N} $

it's not a Dirac delta, it's a Kroenecker (discrete-time) delta.

$$ \begin{align} \frac{1}{N} \sum\limits_{k=0}^{N-1} \left(-W_N^n \right)^k \ & = \ \frac{1}{N} \sum\limits_{k=0}^{N-1} \left(-e^{j 2 \pi n/N} \right)^k \\
& = \ \frac{1}{N} \sum\limits_{k=0}^{N-1} \left(e^{-j\pi}e^{j 2 \pi n/N} \right)^k \\
& = \ \frac{1}{N} \sum\limits_{k=0}^{N-1} \left(e^{j 2 \pi (n-N/2)/N}\right)^k  \\
& = \ \frac{1}{N} \frac{1 - e^{j 2 \pi (n-N/2)}}{1 - e^{j 2 \pi (n-N/2)/N}} \\
& = \ \begin{cases} \frac{1}{N} N \text{  for } n = N/2 + mM \text{  where } m \in \text{integers} \\ \frac{1}{N} 0 \text{  otherwise} \\ \end{cases}
\end{align} $$

so it's not exactly $ \delta[n-N/2] $ but more like discrete kronecker comb.

$$ \frac{1}{N} \sum\limits_{k=0}^{N-1} \left(-e^{j 2 \pi n/N} \right)^k  =  \sum\limits_{m=-\infty}^{+\infty} \delta[n-N/2+mN]$$


Fourier transform of a damped cosine wave with a linear frequency chirp I want to take the Fourier transform of the following transient signal,
$$f(t) = e^{-t/\tau} \cos((\omega_0 + m t)t)$$, where $m$ is some gradient parameter in units of $\rm{Hz}/s$.
I thought this would be quite straight forward -- although most of my approaches have been made using Mathematica -- which struggles to provide anything useful.

I would have assumed the resultant function would have a Lorentzian-like peak profilein a similar way to if one takes the trivial Fourier transform of
$$f(t) = e^{-t/\tau} \cos(\omega_0t)$$.
Does anyone have any ideas on how I can approach this, or, an alternative to Fourier transforming a damped sinusoidal with a linear (or even non linear) frequency, i.e. a chirp.

For Fourier Transform of the LFM chirp portion - you use the Principle of Stationary Phase (POSP). The POSP essentially says the main contribution in the Fourier Integral comes from the portion of where the derivative of the phase is zero - it assumes that the integral of the oscillating components cancel themselves out.

Using the POSP - the Fourier transform of a LFM chirp is another LFM chirp in the frequency domain. The assumption is not very good in cases where you have a low time-bandwidth product.

For the complete signal, you have to options:

Evaluate the Fourier Transform of the magnitude portion and then convolve with the Fourier Transform of the LFM chirp pulse.
Keep the magnitude in while you are doing the POSP evaluation. You have to check on the assumptions the POSP make. If I recall correctly - it is usually assumes a slowly varying amplitude window.

I believe you can find the use of the POSP for LFM pulses in the signal analysis books by Papoulis, and also in "Digital Processing of Synthetic Aperture Radar Data" by Cummings and Wong.


Using Matlab's filter() function for solving difference equation with initial conditions I am trying to solve a difference equation with initial conditions using Matlab's filter function:

$$y[n]=\frac12y[n-1]+x[n],\qquad y[-1]=1\tag{1}$$

For $x[n]=6u[n]$, calculation of $y[n]$ by hand gives

$$y[n]=\left[12-(6-y[-1]/2)\left(\frac12\right)^n\right]u[n]\tag{2}$$

However, the solution produced by filter is different. Could anybody point out my mistake?

The code is shown below:

disp('Solution of a difference equation:')
% The example difference equation considered is
% y(n) = 0.5y(n-1) + x(n)
num_coef = [ 1 ]; % coefficient of x
den_coef = [1 -0.5]; % coefficient of y
n = 0:4; % Considering five samples
x = 6 * ones(1,5) % x(n) = 6u(n)
init_cond = [1];  % y(-1) = 1
y = filter(num_coef, den_coef, x, init_cond) % Solution returned by 'filter'
yM = 12 - (11/2)*(1/2).^n % Solution obtained manually
plot(y);hold on; % PLots for comparison
plot(yM,'r')

The problem is the definition of the initial condition of the function filter(). Note that according to the mathworks reference page, filter is implemented using a direct-form II transposed structure. If you look at the corresponding diagram, you see that the filter states (called $w_k(m)$ in the figure) are defined after multiplication with the respective filter coefficients.

So in your example if you want to initialize filter with some past output value $y[-1]$, your initial state is given by $y[-1]\cdot 0.5$ because that's the value of the input of the delay element.

y = filter(num_coef, den_coef, x, 0.5)

should give you the expected result.

Note that for the general case there's the function filtic, which computes the necessary initial conditions given past output samples (and, if necessary, past input samples). For your example you would use the command

zi = filtic(num_coef, den_coef, 1)

which results in zi = 0.5, as explained above.


What is the difference between a lag compensation and PI control? A lag filter or compensator has the Laplace transform transfer function
$$G_c(s) = \frac{s+z}{s+p},$$
with $-z < -p < 0$. In practice, the effect of lag compensation in feedback control, is to increase the DC gain of the open-loop system by the factor $G_c(0)=z/p$, which improves the ability of the closed-loop system to track low-frequency reference signals. The tradeoff is that the open-loop system incurs a phase shift of $-\pi/2$ at frequencies $z < \omega$.
However, in the limit as $p\rightarrow 0$, we are left with a PI controller, which has the form
$$G_c(s) = K_p + \frac{K_I}{s} = \frac{K_p s + K_I}{s},$$
which has infinite DC gain ($|G_c(0)|$), meaning the closed-loop system now achieves zero steady-state error in response to a step reference. The only tradeoff I can see compared to the lag filter is that there is a phase shift of $-\pi/2$ for $\omega < p$, but I don't see why this is a problem in practice.
Since the steady-state error is improved as the pole moves closer to the origin, why would you ever use a lag filter (with $z > 0$) rather than just place the pole at the origin?

Well, the two systems differ only at low frequencies. In fact, if you define
$$
R_L = \dfrac{\tau_zs+1}{\tau_ps+1},\qquad R_I = \dfrac{\tau_zs+1}{\tau_ps}
$$
you have that $R_I(j\omega)\underset{\omega\to\infty}{\to}R_L(j\omega)$.
Therefore, you can expect that the behavior of the closed-loop system differs only at lower frequencies.
What does the PI do?
Consider a plant $G(s)$, define $L_I(s)=R_I(s)G(s)$ and $L(s)=R_L(s)G(s)$. Define the functions
\begin{align*}
F(s) &= \dfrac{L(s)}{1+L(s)},& S(s)&=\dfrac{1}{1+L(s)}, & D(s)&=\dfrac{G(s)}{1+L(s)}
\end{align*}
$F(s)$ is the transfer function from the reference and the output, $S$ is the transfer function from a disturbance acting at the output of $G(s)$ and the system output and $D(s)$ is the transfer function from a disturbance acting additively on the control input and the system output. As you easily see, at frequency $\omega=0$ you have $F(0j)=1$, $S(0j)=0$, $D(0j)=0$. This implies that at front of any constant unknown references and/or disturbances, you always get zero error, provided that the system is stable.
What a lag net cannot do
With a lag controller you do not have the above robustness property. In fact, even if by using a lag net together with an appropriate feed-forward action, you can still have zero steady state error at front of constant references, that is not a robust solution. The reason is simply that the feed-forward control law will (highly) depend on the plant parameters and therefore, as you put an uncertainty on the plant, you lose the zero steady state error property.
With an integral action instead, as long as the system is stable, you always obtain a zero error given constant references and/or disturbances, no matter how large the parameter uncertainty. The only thing you need to ensure is the closed-loop stability, and note that to stabilize a system you don't really need perfect knowledge of the parameters and in fact typically you can tolerate quite high uncertainties.
So why should I use lag controllers?
Well, if you need robust regulation w.r.t. constant references and/or uncertainties there is no theoretical reason to go for a lag controller. I would definitely go for a PI. However, not all the tasks are about regulating a system at front of constant references. Having a pole in the origin to track a sinusoid is not that awesome, is it? The point is, what do you need your controller to do? Every distinct answer to this question needs a different controller, and you typically chose your lag or PI or other control configuration to shape the frequency response of the closed-loop system. Think about this: how would you design a controller to have zero error given references that are sinusoids at frequency $\omega=3$?


FFT Heart Rate Monitor on Arduino Part of my project is to build heart rate monitor from Nano and a PulseSensor. I am having trouble working with the Fast Fourier Transform part of the code.

This code here is from a library (arduinofft) that can be found in the arduino IDE. When testing the arduino the serial printer returns values of about 51 bpm, which is far to low, compared to other devices of about 70bpm.

What i have done:

changed the sampling frequency to 5.6Hz to get the range of 0 to 2.8Hz range for the hate rate
edited the library to only take the last 25 to 128 samples (only when computing the max peak), to off-set the DC. Previously getting readings of 5 BPM

What i don't understand:

the sampling period, if the samples are taken and put into an array every second what does the period do.
Also if samples are taken every second is this enough to get a good reading on a heart rate shouldn't the sample be maybe of 15-30 seconds?
If this is so, should i put another array at the bottom. To capture the results every second upto 15 or 30, print the average then first in first out  the array of new values?

#include "arduinoFFT.h"

arduinoFFT FFT = arduinoFFT(); /* Create FFT object */

#define CHANNEL A0
const uint16_t samples = 128; //This value MUST ALWAYS be a power of 2
const double samplingFrequency = 5.6; //Hz, must be less than 10000 due to ADC

unsigned int sampling_period_us;
unsigned long microseconds;

double vReal[samples];
double vImag[samples];

#define SCL_INDEX 0x00
#define SCL_TIME 0x01
#define SCL_FREQUENCY 0x02
#define SCL_PLOT 0x03

void setup()
{
  sampling_period_us = round(1000000*(1/samplingFrequency));
  Serial.begin(115200);
  Serial.println("Ready");
}

void loop()
{
  /*SAMPLING*/
  for(int i=0; i<samples; i++)
  {
      microseconds = micros();    //Overflows after around 70 minutes!
      vReal[i] = analogRead(CHANNEL);
      vImag[i] = 0;
      while(micros() < (microseconds + sampling_period_us)){
        }
  }

  FFT.Windowing(vReal, samples, FFT_WIN_TYP_RECTANGLE, FFT_FORWARD);    /* Weigh data */
  FFT.Compute(vReal, vImag, samples, FFT_FORWARD); /* Compute FFT */
  FFT.ComplexToMagnitude(vReal, vImag, samples); /* Compute magnitudes */
  double x = FFT.MajorPeak(vReal, samples, samplingFrequency); /* Compute largest peak*/
  double y = round((60/(1/x)));
  Serial.println(y);  //print bpm
  delay(2000); /* Repeat after delay */
}

Thanks for any help

1) The sampling period ($T_{s}$) is the time between each sample of the observed analog signal that is captured and stored in 'Vreal'

2) The total observation time ($T_{obs}$) is then equal to $N T{s}$ with N being the number of samples you captured

Very roughly, In order to properly monitor your signal, you need $T_s$ to be at the very least 2 times lower than the fastest observed phenomenon.
In other words, you need the sampling frequency ($F_s$) to be 2 times greater than the signal's maximum frequency (see https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem).

Also, you need $T_{obs}$ to be great enough to capture enough heartbeats to the spectral analysis done to determine the heart pulse rate actually makes sense.

What you need to know then, is :

1) How much heartbeats you need to capture so the heartbeat can be properly determined. This will be answered by people working in the field of medicine more than people from signal processing

2) How is a heart beat represented in terms of spectral content, and how large is its spectral bandwidth so you can set up your sampling frequency high enough to be able to capture them without missing critical components


Demodulating an oversampled OFDM signal I'm trying to simulate an OFDM signal that has been captured with a receiver that has a much larger bandwidth than the OFDM signal.  I tried the following MATLAB code to simulate an OFDM signal that's 20% of the captured bandwidth and then demodulate it:

Nsymbol = 100;
M = 16;
Nfft = 128;
CPlen = 12;
% random bits...
in_bits = randi([0 1], log2(M), Nsymbol*Nfft);
% use a 16-QAM symbol constellation
symbols = qammod(in_bits, M, 'InputType', 'bit');
% reshape it ofdmmod
parallel_symbols = reshape(symbols, Nfft, Nsymbol);
baseband = ofdmmod(parallel_symbols,Nfft,CPlen);
% oversample by a factor of 5 and then decimate by a factor of 5 again
oversampled = resample(baseband, 5, 1);
back2baseband = resample(oversampled, 1, 5);
symbols2 = ofdmdemod(back2baseband, Nfft, CPlen);
plot(symbols2(:), '.');

When I run the code, I expect that the symbol constellation should be (roughly) the same as the original 16-QAM constellation, but instead I get this:

As a result, I get a BER of about 1% with no noise in my simulation... clearly I am doing something wrong.  What is the right way to demodulate an oversampled OFDM signal in MATLAB?  I tried using the DigitalUpConverter and DigitalDownConverter objects as well, but that was a total disaster.

EDIT: Thanks to Dan's answer below, I realized that the filters that MATLAB applies when downsampling induced a linear phase shift so I was looking at the wrong sample.  When I just do back2baseband = oversampled(1:5:end) I recover the ideal constellation.

The result of your plot is what you would see when you have timing offset error but no frequency error. Given that it is over-sampled you can confirm this by manually adjusting your decision point a few samples either way. Even better plot out an eye diagram such as what I show below with all your samples for 16 QAM and it should be clear from that if this is indeed the issue.

(With as little as 2 samples per symbol you can completely replicate the eye diagram as in the plots below using the MATLAB "resample" command).

Consider the eye diagram on the right showing the real and imaginary components on the upper and lower graph, and if we sample at time T=0, we should see tiny dots at all the expected locations for 16QAM on an IQ constellation such as the one you show. If we shifted the sampling decision points to the right or left (timing offset), I believe the IQ constellation would look similar to what you show.


clock + phase recovery I am working on the implementation phase and time recovery algorithm and have doubts. Could you please explain to me some question?

PLL for time and phase recovery

I have found I can use PLL not only for a carrier phase correction but also for time. Is there difference in implementation?
Reference: Digital communications a discrete-time approach by Rice

what should be first: phase or time recovery?

A receiver needs to synchronize on those two timings in order to demodulate with best performance:

Transmitter symbol rate
Transmitter carrier

Synchronization state of the art way is to use a PLL. You need two loops processing at the same time. The same PLL architecture can be used for both synchronizations but instantiated separately. That is the case for the filtering part of the loop.
Except:

Each loop has a different timing error detector: One detects symbol rate error and the other carrier phase error. Different algorithms for each.
The feedback correction is applied in different areas in the chain: Symbol rate feedback is applied on a resampler, carrier error correction where the frequency translation is done.

Also important, a PLL is a fine synchronization. Which means that if you have a too large offset it won't be able to synchronize alone. A process ahead in the chain should perform the coarse synchronization.
Timing error detector for each loop should be smartly selected. As they can interfere between themselves. For example Gardner when applicable is able to estimate the symbol rate timing error with a carrier phase offset inside.


discrete cosine transformation (dct-II) own implementation with loops I want to implement the DCT-2, the dctmtx:

$$X[k] = \sum_\limits{n=0}^{N-1} x[n] \,\cos \left(\frac{\pi}{N} \left(n+\tfrac12 \right)k \right) ~~~,~~ \;k=0,1, \dots N-1$$

function [ C ] = dct_mtx( N )
for k = 0:N-1
    for n = 0:N-1
       tmp(k+1,n+1) = cos((0.5+N)*k*pi/N);
    end
end

tmp(1,:) = tmp(1,:) / sqrt(2);
C = tmp;
end

If I call the function C = dct_mtx(3); and check the matrix C'*C = I and is not the identity but almost.

What I have to change in the code to get the Identity in the end?

the identity looks like that:

1.50000000000000         -1.66533453693773e-16    -1.11022302462516e-16
-1.66533453693773e-16    1.50000000000000         -2.77555756156289e-16
-1.11022302462516e-16    -2.77555756156289e-16    1.50000000000000

You are missing a $\sqrt{\frac{2}{N}}$ scaling factor. Change the following line:

tmp(k+1,n+1) = sqrt(2/N)*cos((n+0.5)*k*pi/N);

You can also simply use Matlab's built-in function dctmtx.

C = dctmtx(n)


Generalized translation on graph David I.Shuman in "vertex-frequency analysis on graph" claims that,"we generalize one of the most important signal processing tools – windowed Fourier analysis – to the graph setting and When we apply this transform to a signal with frequency components that vary along a path graph, the resulting spectrogram matches our intuition from classical discrete-time signal processing. Yet, our construction is fully generalized and can be applied to analyze signals on any undirected, connected, weighted graph."

In this paper generalized translation operator that allows us to shift a window around the vertex domain so that it is localized around any given vertex, just as we shift a window along the real line to any center point in the classical windowed Fourier transform for signals on the real line.

generalized translation operator:

$$(T_{i}f)(n) := \sqrt{N}(f*\delta_{i})(n) = \sqrt{N}\sum_{\ell=0}^{N-1}\hat{f}(\lambda_{\ell}) \chi^{*}_{\ell}(i)\chi_{\ell}(n) .$$

See the picture

Question

How do I understand translation on the graph?
How to label irregular graphs with high dimensional data?

Thanks.

To understand either of these, you first have to understand the basic premise behind Graph Signal Processing (GSP) which is to map a signal to a graph and then work with it on the "Graph space". This is possible due to certain similarities of classic DSP concepts and Algebraic Graph Theory.

So, it is easier to start from the second question because, before we start applying GSP, we first need a graph.

  How to label irregular graphs with high dimensional data?

The short answer is that this is still an open problem and currently, there are "signals" that are naturally mapping on graphs and others where the mapping is either arbitrary or in some way constructed.

Signals that naturally map on graphs are usually expressed as weights of the graph's edges through a Weight Matrix that is similar to an Adjacency Matrix.

Typical examples are usually items and some form of similarity between them. For example, suppose that you have a set of $N$ time series $X$ and you evaluate their cross correlation. This will result in a Weight Matrix (let's call it $W$) whose $i^{th}, j^{th}$ element ($W_{i,j}$) is the cross correlation between time series $X_{:,i}$ and $X_{:,j}$. (So, $X$ is an $m \times n$ matrix of $n$ time series signals each being $m$ samples long.)

What does this graph look like? It looks like a Clique. In other words, because we have examined all-to-all cross correlatons, all nodes are considered connected with each other. But, the strength of the connection is expressed by some weight.

So, yes, they are all connected, but some are much more closer than others.

For signals that do not naturally map on graphs, you first have to solve the corresponding graph labeling problem. This is generally done in two ways, either by coming up with a function that maps a signal to some graph or arbitrarily.

In the arbitrary case, you select some graph whose order (the number of nodes) is equal to the number of samples in your signal. That graph's nodes can be arbitrarily connected, it could for example be an entirely random graph where there is equal chance for any two nodes to be connected.

This is what the author of the paper that you link is actually doing. They take an arbitrary graph (a road network) and they map on to it an exponential decay signal. How? Arbitrarily. Does it make sense? No, but it illustrates the point they are trying to make about showing the effect of the operators.

(See page 4:"Note that the definitions of the graph Fourier transform and its inverse [...] depend on the choice of graph Laplacian eigenvectors, which is not necessarily unique. Throughout this paper, we do not
specify how to choose these eigenvectors, but assume they are fixed.  The ideal choice of the eigenvectors in order to optimize the theoretical  analysis conducted here and elsewhere remains an interesting open question;  however, in most applications with extremely large graphs, the explicit computation of a full eigendecomposition is not practical anyhow, and methods  that only utilize the graph Laplacian through sparse matrix-vector multiplication are preferred.")

The other way that you can do the mapping is with an intuitive or model fitting (in the sense of optimisation) way.

So, an intuitive way to map a signal to a graph is to put the samples of some $x[n]$ time series on the nodes of a graph that are simply connected as a "line" (so, something looking like $x[0] \rightarrow x[1] \rightarrow x[2] \rightarrow x[3] \ldots \rightarrow x[n]$ ).

And a constructed way is to use optimisation in order to construct a graph whose connectivity represents SOME aspect of your original signal $x[n]$.

Which brings us to the first question:

  How do I understand translation on the graph?

The short answer is that translation on a graph is equivalent to a re-ordering of the edges that effects a new connectivity pattern on the nodes of the graph. In this way, the nodes appear to have "moved" or translated to a different "position".

So now the question is how do you define "position" and to an extent this question is a bit related to the first one because "position" and how you represent the signal are related.

But, here is a very simple example, just to demonstrate a trivial translation.

Say we have this signal: $x = \left\{ 0,1,2,3,2,1,0,1,2,3,2,1,0 \right\}$ and we map it to the "line" graph $G(V,E)$ we saw earlier that looks like $x[0] \rightarrow x[1] \rightarrow x[2] \rightarrow x[3] \ldots $. In other words, we assign $x[0]$ to $v_0$, $x[1]$ to $v[1]$ and so on and we assume that nodes are connected "sequentially" (and cyclically).

The adjacency matrix (or the weight matrix) of this graph is:

$$A = \begin{vmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
                      1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\end{vmatrix}$$

Notice here that I have connected $x[|V|]$ back to $x[0]$ with that last line.

So, how do you "move" nodes around?

In classic DSP, if you wanted to shift things in time, you did something like: $y[n] = x[n+2], n \in \left\{0 .. |x| \right\}$ and the shift is cyclic here. After this, our shifted sequence $y$ looks like: $y = \left\{ 2,3,2,1,0,1,2,3,2,1,0,0,1 \right\}$

Right, so how could we achieve the same, solely utilising the $A$ to effect the same shifting on our graph signal?

Well, that's easy here because instead of having our initial time series $x$ connected as $x[n] \rightarrow x[n+1]$ it is as if we now connect $x[0] \rightarrow x[2], x[1] \rightarrow x[3], x[2] \rightarrow x[4], \ldots$. So basically, the same $A$ as above, only that now the $1$s appear two places to the right of their current position.

Did you notice how we expressed something that happened in the time domain to something that happens in the graph domain?

The key idea here is that we translated the signal by changing the way the nodes are connected. Translation is basically an operator on the connectivity of the graph.

BUT!

Notice here that we said earlier that we assume a "line" graph as the underlying graph for our signal. We made an arbitrary decision. We could have mapped our signal on the road network of some city (as the authors of the paper that you link have done). Then, how do you define translation on that thing?!?".

This is where the Graph Laplacian and Algebraic Graph Theory come into play.

To cut a long story short, the Graph Laplacian is like the Discrete Fourier Transform for signals in the time domain.

It has been the topic of a lot of research in pure mathematics and its eigenvectors and eigenvalues are supposed to return a lot of information about the graph's connectivity structure (for example, whether it contains cycles or not, what sort of lengths of cycles, whether it is completely connected or not, etc).

So, basically, what the authors are working on in the paper that you linked are a translation operator and a "DFT" equivalent operator on the graph laplacian so that you can "translate" nodes around an arbitrary connected graph (not only one looking like a line, it could have any shape) and decompose and recompose the graph connectivity matrix to elementary components no matter how complex the graph is.

You can see now how the representation of the graph and "translation" are connected. The Laplacian of the graph depends on the values of its adjacency (or weight) matrix (i.e. its structure). You map your signal $x[n]$ on the node set of the graph $V$ and you assume (or construct) the edge set $E$. Therefore, any notions of "translation" or "frequency" now depend on the structure of the adjacency matrix.

Therefore, don't try to understand why Fig.7 in the paper that you link looks the way it looks. First of all, the mapping of the signal on the road network is arbitrary and second, the "translation" depends both on the mapping and the connectivity matrix of the road network. Conceptually, this particular example does not have an immediate connection with reality. But at the same time, conceptually it shows you what translation means over a graph signal and a graph that can have arbitrary connectivity.

Perhaps it is easier to think about GSP in terms of linear algebra because at the end of the day, this is what it is all based on.

If we forget about graphs, adjacencies, nodes, edges, mappings, etc for a minute and focus on the Laplacian: The whole point of GSP is to come up with a new representation for $x[n]$ in the form of a matrix. A new "decomposition" if you like, similar to the way the DFT matrix decomposes a signal or similar to the way Wavelets decompose a signal.

In fact, wavelets are part of the family of "constructed" graphs I am talking about earlier. They are basically a matrix. This matrix could also be expressed as a graph. When it is expressed as a graph, it opens the door to "new" ways of working with signals or "new" ways of working with graphs. For more information on this line of thinking please see this paper or this paper and this paper (for methods of discovering graph representations).

Hope this helps.


Equal contribution criteria during image pyramid construction In their paper on Image Pyramids, Burt and Adelson state the following about pyramid generation.

  An additional constraint is called equal contribution. This stipulates that all nodes at a particular level must contribute the same total weight($=\frac{1}{4}$) to nodes at the next higher level. Let $\hat{w}(0)=a$,$\hat{w}(-1)=\hat{w}(1)=b$,$\hat{w}(-2)=\hat{w}(2)=c$. In this case, equal contribution requires that $a+2c = 2b$.

How is this weight condition derived and why is the weight $\frac{1}{4}$? For a larger context of the question, the screenshot of the relevant section in the paper can be seen below

The following picture from this set of lecture slides explains the equal contribution. The contribution from two nodes in the same layer must be equal. With reference to the figure, the sum of weights on green edges must match the sum of weights on the blue edges.


What happens to the audio of a reconstructed signal if we remove too much frequency? Assuming we have an audio with noise embedded in it, I was wondering what happens theoretically if we remove (filter out) too much frequency when we pass LPF with a cutoff frequency lower than the perfect cutoff frequency needed to reconstruct the signal perfectly? Practically after performing the necessary computations on MATLAB and running the audio, the sound level was somehow lower than that of the perfect reconstructed audio. However, I have no explanation as to why this phenomena occurs.

Theoretically, the error between the original and the reconstructed signal goes up.  Practically, for actual audio, if you're lowpass filtering too much, the audio starts to sound muffled or muddy -- like you're trying to listen to the event through a pillow.
Again in theory, if you've got a bandlimited signal mixed with noise, the best processing for that signal is to low-pass filter it more than the signal is already bandlimited.  There's a lot of theory out there on how to design optimal filters for this.


Phase lag changes cross spectrum I have a same noise passing through two different ARMA filter. The autospectrum of first filter is given by S11=G1*G1*S_in where S_in is input spectrum created by set of random numbers.

Similarly the autospectrum of second filter is given by S22=G2*G2*S_in where S_in is input spectrum created by same set of random numbers as the previous one.

The cross spectrum that i expect is S12=G1*G2*S_in.But the amplitude is reduced is this due to phase shift of the filter or the derivation above is incorrect.

Attached is an image1 explaining the above statements.

Derivation seems to be wrong. I am guessing from derivation that it is power spectral density that you are talking about. Your derivation seems to be wrong, cross density would be $$S_{xy}(x)S_in=G1(x)G2\times(G2^*)(x)S_in$$ or $$S_{yx}(x)S_in=G1\times(G1^*)(x)G2(x)S_in.$$ If you take magnitude spectrum both will be same, phase spectrum will be of 180 shift.


Harris Corners Detector squared derivatives I have trouble understanding something (small) regarding the math behind Harris Corners Detector.
(for example I looked at this explanation - https://aishack.in/tutorials/harris-corner-detector/)

In the equation we get there are squared derivative matrices - $I_x^2, I_y^2$ and also $I_x \times I_y$. I've seen code implementations of this around the web, and it seems to me that these multiplications are element-wise.
My question why? I thought that when you take the square of a matrix it should be matrix multiplication, and not element-wise multiplication...

I'd really appreciate some clarifications regarding this,
Thank you

I think this part in article is starting point of the confusion :

  See how the $I(x+u, y+v)$ changed into a totally different form
  $I(x,y)+ uI_x + vI_y)$?

The author should have write a bit more detailed notation as:

$I(x,y)+ uI_x(x,y) + vI_y(x,y)$

i.e. pointing out that these are element-wise spatial derivatives (like Sobel) not matrix derivatives ( those x and y indicate horizontal and vertical spatial derivative).


Want to do shot boundary detection via SVM, what are some good features? I want to do shot boundary detection via SVM's.

I'm dividing the frame into nxn blocks. Per block I'm finding these features:

shannons entropy
edges (H,V,Diag)
standard deviation

for consecutive frames and corresponding blocks, the above features are compared via differences and L2 norms, I threshold over that, if above threshold then given a 1 else 0. That is, the feature value belonging to that feature for that particular block becomes 1.

This assignment happens per feature for every block. For a frame, the feature vector is the summation of all the feature vectors of its blocks.

Then I take +-N/2 frames around it and concatenate their feature vectors to create a single feature vector (for that frame).

I trained the SVM using this. There are so many normal frames and only a few shot boundary frames. I'm dividing amongst cut,dissolve, fade, other and normal frames.

I'm unable to properly detect the shot boundaries.

Can someone suggest which features I should use and how I should select them for shot boundary detection? I think I may also have not done quite well on the threshold selection, so any ideas how I could do that better?

I'm using TRECVID's 2001/2002 database.

Basically, Shot Boundary detection technique's output is affected by camera motion, object motion and illumination conditions. The feature you choose should be invariant to these issues. The problem in your case may be due to poor training dataset/ imbalanced dataset. threshold can be fixed experimentally and it is tedious.


How do impulse response guitar amp simulators work? I am wondering how impulse response guitar amp simulators/modelers work. I thought it was a matter of convolving a signal of recorded impulse response in time-space with a guitar sample.

I tried to do this with Matlab using conv function. I loaded an example impulse response and recording of a guitar sample and the result was just a badly distorted sound.

I read that using normal convolution will not work and the response of the amp is dynamic and differs based on the input signal amplitude.

I have found this question: Dynamic convolution vs Volterra series, but I do not quite understand how I would apply this to my case and if this even is the right direction.

I am EE student and know a bit about DSP, but we only dealt with basic signal convolutions so I am quite lost on what to try, but I am really curious on how all this works.

When talking about modeling, there are two things that usually get modeled: 1. the guitar amp, and 2. the speaker cabinet. Only the latter is modeled by an impulse response, which means that the cabinet is simply represented by an LTI system and implemented by convolution. This is of course an approximation but it works fairly well. You can find a lot of such measured impulse responses on the internet.

As for the guitar amps, there is no such thing as an "impulse response guitar amp simulator". As mentioned in the comments and in Marcus Müller's answer, guitar amps are much too complicated to be modeled by a simple impulse response. They are non-linear and dynamical, and the corresponding digital models are quite complex. There's a good reason why the big players in amp modeling (Kemper, Fractal Audio Systems, Line6) do not publish their modeling algorithms.

There are two basic approaches to modeling of nonlinear analog circuits: 1. black-box models, which only try to emulate a measured input-output relationship using a set of test signals, and 2. white-box models, which try to simulate the actual analog circuit in all details.
You can read a good summary of the basics in A Review of Digital Techniques for Modeling Vacuum-Tube Guitar Amplifiers.


Why impulse response Why do people calculate impulse response? Generally systems are designed after modelling the possible input signals and calculating the appropriate filters for the desired outputs. Where does impulse response calculation become important?

Having the impulse response allows one to simplify the input of a sampled linear time-invariant system down to one single unit-magnitude input sample.

It is usually fairly simple to build any time domain input from a single unit input (with a linear combination thereof), and thus, given the impulse response for that single input, the output for any time-domain-described input for that system.  Other basis vectors are possible, but may require more or less work (computational effort) for a given type of description of the input.


Zeropadding and amplitude scaling Could somebody tell me what is the DFT scaling of a $N$-length signal which is padded with $M$ zeros ? Is it $1/(M+N)$ or $1/N$, i.e.
$$
Y_k=\frac{1}{N}\sum_{n=0}^{N-1}y(n)e^{\frac{2\pi nk}{(M+N)}} \quad \text{or}\quad
Y_k=\frac{1}{(M+N)}\sum_{n=0}^{N-1}y(n)e^{\frac{2\pi nk}{(M+N)}}?
$$

The conventional definition of the DFT for a length $N$ signal (without zero-padding) is

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{1}$$

So there is no scaling involved. Scaling is applied to the inverse DFT:

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{2}$$

If you pad a length $N$ signal with $M$ zeros, the resulting DFT length is increased to $N+M$, and the DFT becomes

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/(N+M)}\tag{3}$$

The corresponding inverse DFT is

$$x[n]=\frac{1}{N+M}\sum_{k=0}^{N+M-1}X[k]e^{j2\pi nk/(N+M)}\tag{4}$$

EDIT:

In reaction to your comments, it is indeed necessary to sum over all indices $k$ up to $N+M-1$ in the sum of Eq. $(4)$. A simple example will illustrate this. Take a sinusoidal signal with an integer number of periods inside the DFT length:

N = 16;
n = 0:N-1;
x = sin(pi/4*n);
X = fft(x);

Its DFT X consists of two impulses. However, the DFT of the zero padded signal shows the underlying sinc-like behavior due to the rectangular window. Note that the DFT is just a sampled version of the DTFT, which has a continuous frequency variable. Zero padding just increases the sample density in the frequency domain:

xz = [x,zeros(1,N)];
Xz = fft(xz);
subplot(1,2,1),stem(abs(X))
subplot(1,2,2),stem(abs(Xz))


image resolution in smartphones and laptop If I used a smartphone with high pixel per inch (ppi) value for capturing an image with certain dimensions, what happened if this image is displayed on a laptop (larger screen) with:

the same ppi value
lower ppi value
higher ppi value

Do the dimensions of the image change? or still the same?
Is interpolation required when displaying an image on a mobile phone different from the mobile phone that captured it, if both mobile phones have different ppi values?

Based on your comment:

  I mean that the image is captured using say, for example, the 5
  megapixels camera of the smartphone, and displayed on its screen
  having a certain dpi value. What happens if the same 5 megapixels
  image is displayed on another smartphone having different screen size
  and different dpi?

If you display an image (taken by whichever device, with whatever megapixel camera) on two different displays with two different DPI or PPIs, then it will seem smaller on the high PPI screen and larger on the low DPI screen.

This is very easy to see. Lets say the image has 2000V x 3000H pixels. And display this image on a screen with 200 PPI. Then it will take 15 inches horizontally and 10 inches vertically. On the other hand, if you display it on a second device with 100 PPI, then it will take 30 inches horizontally and 20 inches vertically, hence enlarged by 4 times. Lets conclude :

  Everything else being the same except the screen PPIs; a given image displays larger on a low PPI screen
  and smaller on the high PPI.


Are ACF plots sufficient to detect white noise assuming probable existence of non linear associations? ACF plots use autocorrelations tests to determine associations in time series. Autocorrelations, however, only measure linear associations. To my understanding, certain non-linear associations can exhibit zero autocorrelations. Therefore, can there be a scenario where an ACF plot that is considered to reflect white noise be actually representative of non-linear associations? In other words, are ACF plots sufficient for determining whether a signal is a white noise?

In other words, are ACF plots sufficient for determining whether a signal is a white noise?

A plot is never sufficient. An impulse-shaped ACF is. (A plot can only ever be an estimate based on observation; to find the true ACF, you'd have to observe the process for infinity, which you can't.) You can only observe that the plot approximates the ideal ACF sufficiently much (to your definition of "sufficient").
Anyway, assuming you mean "ACF being a dirac delta", yeah, that fulfills the definition of white. That is sufficient, because it literally is the definition.

can there be a scenario where an ACF plot that is considered to reflect white noise be actually representative of non-linear associations

No, because white noise is defined to be any random process with a delta dirac ACF.

certain non-linear associations can exhibit zero autocorrelations

Yep! "associations" is a bit of very vague term, though. I'd prefer if you said:

Are there discrete random processes (i.e. random sequences) where elements aren't independent from each other, but still be uncorrelated?

Independence is a strictly defined term: It means that the joint density functions of multiple points in time is the product of the individual densities at these times. I.e.: Let $X(t, \xi)\in \mathbb C, \, t \in \mathbb Z$ be a process at times $t$ (with different realisations $\xi$).
Then, the ACF is defined as
$$\phi_{XX}(t_1,t_2) = \mathbb E\left(
X(t_1)X^*(t_2)
\right)$$
and a process is white iff (assuming wide-sense stationarity, so that $\phi(t, t+\tau) =: \phi(\tau) \forall t$)
$$\phi_{XX}(t_1,t_2)=\phi_{XX}(\tau)=\begin{cases}
> 0 & t_1=t_2\iff \tau=0\\
= 0 & t_1\ne t_2\iff \tau \ne 0
\end{cases},$$
from which it's easy to see that different elements from an observed series are uncorrelated.
However, as you probably know, correlation is not the same as independence.
Counterexample:
$$Y(t) = \begin{cases}
U(t) & t\ge 0\\
S\cdot Y(-t) & t < 0
\end{cases},$$
with $U(t)\sim\mathcal U([-1,1])$ (independently uniformly drawn from the $[-1,1]$ interval), and $S\sim\mathcal U(\{-1,1\})$ (a randomly independently drawn sign).
Obviously, different elements from the sequence are not independent. $Y(-5)$ is either $Y(5)$ or $-Y(5)$. Other values are not allowed.
Checking the ACF, you'll see
\begin{align}
\phi_{XX}(t_1,t_2) &= \mathbb E\left(
X(t_1)X^*(t_2)
\right)\\
&=\begin{cases}
\mathbb E\left(
U^2
\right) =\frac13&\text{for }t_1=t_2&\text{(Variance of uniform continous RV)}\\
0 &\text{for } t_1\ne -t_2 &\text{(the $U(t)$ are uncorr.)}\\
\mathbb E\left(
S\cdot U^2
\right)
 &\text{for } -t_1=t_2
\end{cases}\\
&=\begin{cases}
\frac13&\text{for }t_1=t_2\\
0 &\text{for } t_1\ne -t_2\\
\mathbb E(S)\mathbb E\left(
U^2
\right)
 &\text{for } -t_1=t_2 & \text{(Since $S(t)$,$U(t)$ indep.)}
\end{cases}\\
&=\begin{cases}
\frac13&\text{for }t_1=t_2\\
0 &\text{for } t_1\ne -t_2\\
0\cdot\frac13
 &\text{for } -t_1=t_2 \\
\end{cases}\\
&=\begin{cases}
\frac13&\text{for }t_1=t_2\\
0 &\text{for } t_1\ne t_2
\end{cases},
\end{align}
which seems to say "this process is white".
BUT: this process is not Wide-Sense Stationary, so applying the term "white" to it doesn't "work".  So, you need to be very careful to what you apply the label "white", to begin with.


Weight vector of an exponential moving average? With weight vector I mean the vector with weights that you have to multiply the observations in the window that slides over your data with so if you add those products together it returns the value of the EMA on the right side of the window.

For a linear weighted moving average the formula for finding the weight vector is: (1:n)/sum(1:n) (in R code). This series of length n adds up to 1. For n=10 it will be

  0.01818182
  0.03636364
  0.05454545
  0.07272727
  0.09090909
  0.10909091
  0.12727273
  0.14545455
  0.16363636
  0.18181818

the numbers 1 to 10 / 55, with 55 the sum of the numbers 1 to 10.

How do you calculate the weight vector for an exponential moving average (EMA) of length n?

Wikipedia  tells me:

  if n is the length of the window, then alpha<-2/(n+1) and i<-1:n so EmaWeightVector<-((alpha*(1-alpha)^(1-i)))

Is this correct?

Even though the EMA is not really confined to a window with a start and an end, shouldn't the weights add up to 1 just like with the LWMA?

ok the solution was: (in R)

emaweights<-function(m)
{
  alpha<-2/(m+1)
  i<-1:m
  sm<-sum((alpha*(1-alpha)^(1-i)))
  return(((alpha*(1-alpha)^(1-i)))/sm)
}

Here the weights do add up to 1.


Difference between two Sine wave equations I'm confused as what is the difference between these two sine wave equations

x(t) = Sin(2*pi*F*t)

and

x(t) = Sin(2*pi*F*nTs)

x(t) = Sin(2*pi*F*t) is supposed to be used with analog signal, here t is supposed to me the samples when speaking of matlab.

A = 4;       %Amplitude
Freq = 100;
Time = 1/Freq; %0.01
Fs = 1000;      %Sampling Frequency
Ts = 1/Fs;   %Sampling Rate
t = 0:Ts:(Time);
x=A*sin(2*pi*Freq*t);

The output,

Similarly if I use the stem command than a plot command I will get the sampled values instead.

If I'm getting the sampled result from the first equation than what is the purpose to replace t by nTs in the equation and using the second equation ?

I couldn't be able to demonstrate the difference and visualize it in Matlab. Kindly guide me.

The first equation:

x(t) = Sin(2*pi*F*t)

is for a continuous ("analogue") signal, where t may take any real value.

The second equation:

x(t) = Sin(2*pi*F*n*Ts)

(I assume that this is the correct version - it looks like you have a typo in the question ?)

is for a discrete ("sampled" or "digital") signal, where n is the sample number (integer) and Ts is the sample interval (inverse of sample rate). It has values only at the discrete sampling points where t = n * Ts.

In your MATLAB code you are effectively sampling a continuous function, so although you start off with the first (continuous) equation you end up with the second equation because you are only evaluating (sampling) the continuous function at a set of discrete points.

Note also that as Ts -> 0 the discrete version tends towards the continuous version.


Modulation GMSK with I/Q I have a doubt about something. GMSK Modulated signal is equal to I+Q (with I = cos and Q = sin) or I+i.*Q ?

Thank you.

Best regards.

As far as I know, FSK-type modulations do not use quadrature. My answer below is about quadrature, not about a specific modulation.

There are two ways to look at quadrature (I and Q) modulation: baseband, and bandpass.

In baseband, the signal is complex. You have two independent signals, each modulated using PAM, $s_I(t)$ and $s_Q(t)$. The quadrature signal is $$s_{BB}(t)=s_I(t)+js_Q(t).$$ Naturally, this signal can't physically exist. To transmit it, you need to upconvert $s_{BB}(t)$ to some carrier frequency $f_c(t)$. It becomes a bandpass signal: $$\begin{align} s(t)&=\Re(s_{BB}(t)e^{j2\pi f_ct})\\ &=s_I(t)\cos(2\pi f_ct)-s_Q(t)\sin(2\pi f_ct). \end{align}$$

Both signals carry the same information and are equivalent in many ways. $s_{BB}(t)$ is useful to simulate and analyze the system, because it does not have a carrier frequency. $s(t)$ is the signal that is actually, physically transmitted.


How does Fourier Transform manipulates image technically? I need a little bit more explanation on technical details regarding Fourier Transform in Image Processing, mainly on the frequency aspect.

I see that we can use Discrete Fourier Transform to change an image into a signal. From spatial domain to frequency. But from this definition, suppose I want to manipulate, let's say, for example, overall intensity of an image, once my image is in frequency domain, how can I know which frequency is what in the spatial domain? How could I know which frequency belongs to what part of an image?

If there is anything seems incorrect, please correct me.

An answer to another question here provides a great explanation (with pictures) of what a single frequency looks like in the spatial domain.

Note that a single sample in the frequency domain affects all pixels in the spatial domain.  That frequency produces a sinusoid with one frequency in the x dimension and another in the y dimension.

If you want to change the overall intensity of the image, you can do so by manipulating the DC value in the frequency domain.  This is the value at zero frequency in x and y (upper left sample in the frequency domain).  The DC value represents a wave with zero frequency in the spatial domain, so it is just a constant offset over the entire spatial image, which is basically the overall intensity.

If a single frequency affects the entire spatial image, how do you manipulate the frequency samples to do meaningful things in the spatial domain?  Well, one way to think about it is that hard edges and quick intensity changes in the spatial domain correspond with higher frequencies while smooth changes in intensity in the spatial domain correspond with low frequencies.

So, if you want to emphasize edges in the image (sharpen), increase the magnitudes of higher frequencies.  Similarly, if you want to smooth the edges (blur), decrease the magnitudes of higher frequencies.


scipy.signal.decimate vs array operations The following program generates a signal with a bit of noise according to some bits.  I then try to decode the signal using convolve and decimate.  2 problems I notice:

The whole this gets thrown out of wack depending on the 'type' parameters that I pass to convolve, is there a better way of lining up the signal?
The decimate function seems to mess things up - everything works better if just use array indices to access the data I want.  Am I doing something wrong with decimate?

Here is a picture of the output:

import matplotlib.pyplot as plt
from math import sin, pi
from numpy import ones, zeros
from numpy.random import random_sample as sample
from scipy.signal import decimate, convolve

N = 100
amplitudes = [1, 0, 0, 1, 0, 1, 1, 0, 1, 0]
pulse = ones(N)

#generate the transmit signal
tx_pulse = pulse + [0.2*sin(n) for n in range(1,N+1)]
signal = zeros(N * len(amplitudes))
for i in range(0, len(amplitudes)):
    for j in range(0, N):
        signal[i*N+j] = amplitudes[i] * tx_pulse[j]
signal = signal + 0.01*sample(len(signal)) #add noise

#now try to get the data back...
signal2 = convolve(signal, pulse/N, mode='valid')
signal3 = signal2 > 0.5 #bit decision

#dec = decimate(signal3, N, ftype='fir') #THIS MESSES UP THE THE FIRST SYMBOL
dec = signal3[0:N*len(amplitudes):N]

plot4 = zeros(N*len(amplitudes))
for i in range(0, len(amplitudes)):
    plot4[i*N] = dec[i]

print dec

plt.plot(signal)
plt.plot(signal2)
plt.plot(signal3)
plt.plot(plot4, 'ro')
plt.show()

You should use the 'full' parameter for the convolution. Then, naturally with the convolution, your peak that you use for decimation occurs at the end of each symbol. At this point, the transmit and receive filter line up and you get the maximum overlap and hence best SNR.
The decimate function applies a low-pass filter before the decimation. This is due to the sampling theorem: When you reduce the sampling rate, you need to remove higher frequencies to prevent aliasing. However, in your case, you have already convolved with your filter, so the simple version of taking only every n-th sample is perfectly correct.


Finding the frequency response of a filter defined as a Z-domain transfer function I am software engineer studying DSP on my own. I came across this question in a graduate textbook and am not sure how to proceed.

Determine the frequency response for the filter:
H(z) = 1.5 / (1 + 3z^-1)

Clearly, this is a z-domain transfer function. How do I get the frequency response? I am thinking of converting it into a time-domain difference equation, calculating about 128 points, running it through MATLAB's FFT, plot it, and voila, the final curve would be the answer. Is that a right approach? Is there a more analytic/math-y answer?

The analytic way is to substitute the variable $z$ by $e^{j\omega}$ to get the frequency response $H(\omega)$ (with $\omega = \frac{2 \pi f}{F_s}$) - that is to say, the frequency response is the $z$ transform evaluated on the unit circle.

Note that matlab has a built-in function for plotting the frequency response straight from filter coefficients (freqz), whose source code can be seen. Have a look! The implementation does not even involve solving the difference equation or filtering an impulse - you can simply divide the FFT of the numerator sequence by the FFT of the denominator sequence - padded up to the desired number of frequency points.


Wiener filter to estimate constant signal I want to compare the performance of a Wiener Filter and the Kalman filer to estimate the value of a constant $d$ using mesurements corrupted by a white noise. That is, my measurements are of the form
$$x(n) = d + v(n)$$
where $v(n)$ have a normal distribution with mean $0$ and known variance $\sigma^2$.

Using the Kalman filter, I could put it on State Space form
$$d(n+1) = d(n)$$
$$x(n) = d(n) + v(n)$$
and solve the problem. But I am having difficults to set the problem so I can solve with the Wiener filter. My desided signal is the constant signal $x$. My filter input is the measurements $x$. But I have doubts if I am doing something wrong, because my estimated signal is not a really good estimative of the desired response. I used the following matlab code:

function test()
    n = 0:511;
    d = 10 * ones(1,512);
    v = 0.5*randn(1,512);
    x = d + v;
    w = WienerFIRFilter(x, d, 12);
    y = filter(w', 1, x);

    plot(x)
    hold on
    plot(y, 'r')
end

The function WienerFIRFilter is defined as following

function w=WienerFIRFilter(u,d,M)
    aux = xcorr(d,u,'biased');
    p = aux(1,(length(aux)+1)/2:((length(aux)+1)/2)+M-1);
    [U, R] = corrmtx(u,M-1);
    w=inv(R)*p';
end

Am I doing something wrong?

You can use this code

 function w=WienerFIRFilter(input,desired,M)
     auto_corr=xcorr(input,M,'unbiased');
     r=auto_corr(M+1:end); % positive lags only
     R=toeplitz(r); %correlation matrix
     cross_corr=xcorr(input,desired,M,'unbiased');
     p=(cross_corr(M+1:end));
     w=inv(R)*p'
 end


Question about path loss exponent estimation from measurment I have a simple question that confuses me.

Assume that we have the following results of a measurement campaign, distance (km) and the corresponding path loss $L_{e}$ (dB). I want to estimate the path loss exponent of the simplified path loss model

$$L_{m,\text{dB}} = K_\text{dB} - 10 n \log_{10}\left(\frac{d}{d_0}\right)\text ,$$

where $K$ is a constant that depends on antenna characteristics and average channel attenuation at reference distance $d_0$ and is given by $K_\text{dB} = 20 \log_{10}\left(\frac{\lambda}{4 \pi d_0}\right)$ or $K_\text{dB} = 20 \log_{10}\left(\frac{\lambda \sqrt{G_t G_r}}{4 \pi d_0}\right)$, where $G_t$ and $G_r$ are the transmit and receive antenna gains, respectively, and $\lambda$ is the wavelength.

To estimate the path loss exponent $n$, I am finding $n$ that minimizes the MMSE error for the dB power measurement

$$\begin{align}
\text{MMSE}(n) &= \sum_i^N \left[L_e(d_i) - L_m(d_i)\right]^2\\
 &= \sum_i^N{ \left[L_e(d_i) - K + 10 n \log_{10}\left(\frac{d_i}{d_0}\right)\right]^2}\text ,
\end{align}$$

where $N$ is the total number of measurements.

It seems the exact value of $n$ depends on the reference distance $d_0$ and the antenna gains $G_t$ and $G_r$, right?
In other words, for different $d_0, G_t, G_r$, I will have different value for $n$.  Is this in conflict with the fact that on a log-log scale, the path loss is represented by a straight line with a slope of $10 n$ per decade?

The path loss exponent does not depend on the reference distance. There are two mistakes in your path loss formula. First, path loss does not include the gains of transmit or receive antennas. Second, the path loss exponent must also be part of the expression for the path loss at the reference distance.

According to the Friis transmission equation  the ratio between the received power $P_r$ and the transmitted power $P_t$ is given by

$$\frac{P_r}{P_t}=G_tG_r\left(\frac{\lambda}{4\pi d}\right)^2\tag{1}$$

where $G_t$ and $G_r$ are the gains of the transmitter and receiver antennas, respectively, $\lambda$ is the wavelength, and $d$ is the distance between the transmit and receive antennas. The squared term in parentheses (or its inverse) is the free-space path loss. Note that it does not include the antenna gains. In order to account for other environments than free-space, a path loss exponent $n$ is introduced:

$$\frac{P_r}{P_t}=G_tG_r\left(\frac{\lambda}{4\pi}\right)^2\frac{1}{d^n}\tag{2}$$

In free space we have $n=2$; in lossy environments, $n$ can be greater than $2$.

If we introduce a reference distance $d_0$, Eq. $(2)$ becomes

$$\frac{P_r}{P_t}=G_tG_r\left(\frac{\lambda}{4\pi}\right)^2\frac{1}{d_0^n}\left(\frac{d_0}{d}\right)^n\tag{3}$$

In dB we get

$$\begin{align}10\log_{10}\left(\frac{P_r}{P_t}\right)&=10\log_{10}(G_t)+10\log_{10}(G_r)+\ldots\\&\ldots +20\log_{10}\left(\frac{\lambda}{4\pi}\right)-10n\log_{10}(d_0)-10n\log_{10}\left(\frac{d}{d_0}\right)\end{align}\tag{4}$$

where

$$PL_0=20\log_{10}\left(\frac{\lambda}{4\pi}\right)-10n\log_{10}(d_0)\tag{5}$$

is the path loss at the reference distance $d_0$. Note that in $(4)$ the reference distance $d_0$ will cancel out, so the estimate of the path loss exponent $n$ will not depend on $d_0$. On the other hand, the estimate of $n$ will have a (weak) dependence on the antenna gains. This is however not in conflict with the fact that the path loss as a function of distance is a straight line on a log-log plot.

Note that in practice, the path loss at reference distance $d_0$ is often just set to the free-space path loss, i.e using $n=2$ in Eq. $(5)$ $[1]$. In that case the reference distance $d_0$ will no longer cancel out in $(4)$ and the estimate of the path loss exponent $n$ will depend on the choice of $d_0$. Alternatively, $PL_0$ can also be measured at distance $d_0$, or its can be optimized together with $n$ to obtain an optimal fit to measured data $[1]$.

[1] Wireless Communications, A. Goldsmith, 2005 Cambridge University Press, p. 41.


Recognizing hand-written archaeological signs I have 75 images of handwritten signs from which I extracted 7 Hu moments and solidity features. How can I find similarities among them to train a classifier and predict the value? I thought SVM was a good choice, but I don't have a target vector (what do I put? I do not know differences in signs that I can say there are, say, three labels, such as "circle", "triangle" or "square").

Is cluster analysis more appropriate? If yes, What is best method?
I am using the scikit-learn module in python.

http://scikit-learn.org/stable/modules/svm.html

http://scikit-learn.org/stable/modules/clustering.html

It is a bit hard to understand what you are trying to do. What are these signs? The one you posted looks like a wheel.  Are there meaningful categories that you can name? If so, then this is a supervised learning (classification) problem, and you should use a classification algorithm such as SVM.

If there are no clear labels, but you want to group together similar-looking signs, then this is an unsupervised learning problem, and you should use a clustering algorithm.

From your description it sounds more like a clustering problem.


Help generating coefficients for basic FIR LPF I'm trying to generate coefficients for a FIR low pass filter to be used in an FPGA.  I'm using python (scipy.signal) to attempt to do this, but am having trouble getting coefficients in a usable form.

This is the function I'm using to generate coefficients:

b = signal.firwin(21, 0.01, window='blackman');

The coefficients it returns are:

[ -1.62779573e-18   1.08170091e-03   4.74493669e-03   1.19927580e-02
   2.37996986e-02   4.03774280e-02   6.06306778e-02   8.20599398e-02
   1.01201426e-01   1.14487974e-01   1.19246920e-01   1.14487974e-01
   1.01201426e-01   8.20599398e-02   6.06306778e-02   4.03774280e-02
   2.37996986e-02   1.19927580e-02   4.74493669e-03   1.08170091e-03
  -1.62779573e-18]

These are unusable to me when it comes to FPGA implementation.  How do you scale these or is there a better technique that produces either integer coefficients or coefficients that could be scaled to fixed point in hardware?

These even aren't great coefficients for any software implementation that supports floating point – you could easily omit the last and the first coefficient without changing the filter significantly; the rest of the coefficients seem relatively well-conditioned.

So, basically, what you do is that you pick a fixed point bit width $B$ that suits the signal you're dealing with and the precision you need, and multiply each coefficient with $k=\frac{2^{N}}{\max\limits_{x\in b}|x|}$ (or $2^{N-1}$, as we're probably talking about signed numbers); you then usually just round, and try your filter with pseudorandom white noise, as you'll be inferring some degree of quantization error. Alternatively, you analytically calculate the frequency response in $\mathcal Z$-domain (freqz does that for you) of the quantized filter. Of course, you'll see an additional gain of $k$, but DSP engineers typically just "carry through" such factors (as most operations are linear), and interpret the results accordingly at the very end.

Now, as said, your coefficients (aside from the edge ones) look relatively well-conditioned, and I daresay that if you for example decided that 24 bit are enough precision for your coefficients, then that quantization error might be negligible.

If, however, the resulting frequency response is suboptimal, you can increase bit depth of the coefficients in bit steps¹, or by selectively increasing or decreasing symmetrical pairs of taps. For example, if you find that your frequency response is distorted to have too much gain at low frequencies, it might be a good idea to alternatingly in- and decrease adjacent coefficients by one (be slow with that, and always make sure the errors you're fixing even matter at all – who cares if your passband ripple is 0.2 dB instead of the floating point 0.18 dB, if all you needed was 0.4 dB?).

Without loss of generality, yes, there's also methods to generate integer taps directly (they shouldn't be much better for the window method you're using than rounding afterwards, but might be better for iterative approaches), but I'm not versed in these, or would even know of any toolbox that does them.

¹: as FPGAs really do not care whether you're using 23, 24, or 25 bits – unlike CPUs, they are not inherently limited to specific bitwidths; now, this is a bit of a lie, as for multiplication you might be using "DSP slices", which do have finite bit width and are fixed, not implemented on demand. However, it's thoroughly possible that your synthesizer might decide that for "fixed factor" multiplications, a straightforward combinatorial multiplier is appropriate. Maybe don't worry about these details right now.


Normalizing the magnitude of Goertzel Filter c# I'm trying to detect some frequencies of 18Khz with sample rate of 44100 [samples/sec].

When I pass a buffer(~4500 samples) to the algorithm it returns a magnitude in a very strange range. Sometimes i get 0.1, or 7 or 250(!) or 154 or 25000 - very strange. I can't decide on a threshold boundary for the filter. Is there any way noramlizing that result? (My result is the magnitude of the complex calculation of real and image)

This is my Goertzel code:

public float goertzel_mag(float[] data,int numSamples, int TARGET_FREQUENCY, int SAMPLING_RATE,int start,int end)
    {
        int k, i;
        float floatnumSamples;
        float omega, sine, cosine, coeff, q0, q1, q2, magnitude, real, imag;

        float scalingFactor = (float)(numSamples / 2.0);

        floatnumSamples = (float)numSamples;
        k = (int)(0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
        omega = (float)(2.0 * Math.PI * k) / floatnumSamples;
        sine = (float)Math.Sin(omega);
        cosine = (float)Math.Cos(omega);
        coeff = (float)2.0 * cosine;
        q0 = 0;
        q1 = 0;
        q2 = 0;

        for (i = start; i < end; i++)
        {
            q0 = coeff * q1 - q2 + data[i];
            q2 = q1;
            q1 = q0;
        }

        // calculate the real and imaginary results
        // scaling appropriately
        real = (q1 - q2 * cosine) / scalingFactor;
        imag = (q2 * sine) / scalingFactor;

        magnitude = (float)Math.Sqrt(real * real + imag * imag);
        return magnitude;
    }

Your Goertzel code looks okay, the standard fare.

Regarding normalization and thresholds, if you're detecting 18kHz, you should calibrate your signal at different levels and establish a threshold from that data

Start with full-scale (1.0, 0dBFS) 18kHz tones, as those should give you the highest magnitude over the course of your buffer (4500 samples) and then
Repeat the Goertzel filter mag measurements on tones all the way down in level, as much as you care but I would measure at least 6dB.
Finally, run your Goertzel filter over non-tonal material, or buffers of audio that don't contain what you're looking for. This will give you an idea of what the filter will output in the non-tone case. In my application, environmental audio/speech can have incidental tones in the range I'm looking for (400-600Hz) and I want to make sure to understand what Goertzel magnitudes I would get versus Goertzel magnitudes in the presence of a real tone.


the $L^2$-norm of a signal is also applied as its energy! I am  a  newcomer in signal processing. I saw that the $L^2$-norm of a signal is also applied as its energy! How is this concept illustrated for those ones who are working in pure math.

Yes, the square of the $L_2$ norm of a signal is also by definition its energy $\mathcal{E}_x$.
There's nothing surprising, unbelievable, or mysterious in that though?
The concept of signal energy :
$$ \mathcal{E}_x = \int_{-\infty}^{ \infty } x(t)^2 dt\tag{1} $$
is fundamentally based on the concept of energy (or work) in physics as the Kinetic Energy of a particle with mass $m$ and velocity $v$ given by
$$ K = \frac{1}{2} m v^2 \tag{2}$$
There is also the concept of power defined as the rate of work $W(t)$ (or energy).
$$ p(t) = \frac{dW(t)}{dt}  \tag{3} $$
The relation between instantaneous power $p(t)$ and the total energy is
$$ \mathcal{E} = \int_{-\infty}^{\infty} p(t) dt \tag{4} $$
Electrical engineers ignore the mechanical roots, and rely on an electrical analog of energy as heat loss in an Ohmic resistor defined to be:
$$ \mathcal{E} = \int_{-\infty}^{\infty} p(t) dt \tag{5} $$
Where $p(t)$ is the instantaneous electric power associated with a current $i(t)$ passing through a linear time-invariant resistor $R$ , and is given by :
$$ p(t) = R  \cdot i^2(t) \tag{6} $$
( $p(t) = v^2(t)/R $ is also an equivalent expression, based on Ohm's law $v(t) = R i(t)$)
Then the energy of the current signal passing though a linear time invariant system (denoted by simple Ohmic resistor $R$) is given by :
$$ \mathcal{E} = \int_{-\infty}^{\infty} R \cdot i^2(t) dt \tag{7}$$
Now at this point, just ignoring the resistor $R$ (or setting it to $R=1$), and replacing the current variable with a general unitless $x(t)$, we arrive at the mathematical definition of signal energy of as:
$$ \mathcal{E} = \int_{-\infty}^{\infty} x^2(t) dt \tag{8}$$
That being clear, in a parallel course, the study of normed linear Hilbert spaces also consider mathematical p-th Euclidean norm of a complex valued vector as :
$$ L_p = \left( \int_{-\infty}^{\infty} |x(t)|^p dt \right)^{1/p} \tag{9}$$
And you can see that the square of the case $p=2$ corresponds to the signal energy as defined in Eq.(8).
All of these can also be transferred to discrete-time domain.


How many 8x8 JPEGs exist with distinct DCT given a quality level? For a given quality of image (or quality of compression), Multiple 8x8 JPEG images can have the same DCT (Discrete cosine transform) matrix, for a given quality of JPEG compression. I want to find the number of those 8 by 8 JPEG files which have different DCT codes.

I write the question in a different way:

  Assume a set $A$ of images, which is initially empty. We consider all possible DCT's for 8 by 8 Jpegs one by one, and convert each one to an image. If the image already exists in the set $A$, we throw it away, otherwise, we put it in set $A$. After we are done with all possible DCT's, how many images are in the set $A$? (what is the cardinality of $A$?)

Wikipedia has the formula for the DCT conversion. But I cannot realize from that how to find what I need. I am wondering how to solve it or if that is already solved in the literature. Any idea?

Well, $8\times 8$ is simply 64 individual values.

Now, if you had a quantization matrix $Q$ (dimensions $8\times 8$) with all $1$s, you'd have 8 bit to spend on every entry in your DFT, so that's

$$(256)^{8 \times 8}=({2^8})^{(8 \times 8)} = 2^{512}$$

possible different DFTs.

Now, typically, that's not happening, since quantization matrices restrict the possible values that can end up being non-zero.

The number of possible values is relative easy to figure out: You take the maximum value for each entry, 255, and divide it by the quantization matrix value, round that down; that leads you to the number of possible combinations being

$$ N = \prod_{i = 1, j = 1}^{8,8} \left\lfloor\frac{255}{Q_{i,j}}\right\rfloor$$


Spectral Property of Bandlimited Quantized Signal recently I meet a problem.
I generate a random Gaussian signal $x$ with mean $0$ in MATLAB. And use a $4$th lowpass Butterworth filter to acquire $x'$. Then I use the binary rounding quantization method to get the quantized signal $x''$. And I find that the Fourier spectrum (FFT) of $x''$ contains some high-frequency information. That means some energy leaks to high-frequency regions after quantization.
Also, I find the energy in high-frequency areas is related to the bandwidth of the filtered signal $x'$. When the bandwidth of the filtered signal $x'$ is relatively low,  the energy of the high-frequency areas of the quantized signal $x''$ is also lower. I find that some historical references about the assumption of the whiteness for quantization noise. But it seems that this is not the case.
Do you all meet this question? Could anyone give me references for this problem? Thanks in advance!

I guess this probably just a mistake in your analysis code.
Quantization noise is white and for a noise signal it's uncorrelated to the original signal, so spectrum of the original noise doesn't matter.
I did repeat your steps and saw exactly what I expected: The quantization noise is white and the spectrum of the quantized signal follows the original signal at low frequencies and is white at high frequencies (since it's dominate by the quantization noise).
Potential errors in your code

Clipping the signal
Improper spectral estimation, i.e. wrong  choice of window and/or overlap. Not windowing or properly during spectral analysis. Either window manually or just use PWELCH instead.

Here is a MATLAB code snippet that shows this
%% script to generate quantization noise
nFFT = 16384;
fs = 48000; % sample rate
fc = 4000;  % cut off for low pass filter

% create noise, filter, and normalize
x0 = randn(10*nFFT, 1); % 10 frames or so
[b,a] = butter(4,fc*2/fs);
x0 = filter(b,a,x0);
x0 = x0./max(abs(x0));

% quantize to 8 bit
xq = round(2.^7*x0)*2^(-7);

% calculate quantization noise
noiseSignal = xq-x0;
% power spectrum
yAll = [x0 xq noiseSignal];
psd = pwelch(yAll,hamming(nx));
freqAxis = fs*(0:nFFT/2)'/nx;
% graph
clf; figure;
plot(freqAxis,10*log10(psd));
title('Power Spectra');
ylabel('Level in dB');
xlabel('Frequency in Hz');
legend('Original','Quantized','Noise');
grid ('on');


I am not sure why the first 3 coefficients in python are always different than that of MATLAB.

It's just a scaling. If you call Matlab's zp2sos without the second output argument, you will get the same three values as Python. If you need Python to match Matlab, use the first value of the sos as your "k" and divide the first three values in the sos by k (in your example this yields 1 2 1).

Output ordering for 'p' is different than MATLAB but all values account for. Not sure if it matters?

Depending on your application it actually does matter. For rather complicated reasons, Python's is better and will in most cases generate less numerical noise  and often significantly so.


NxM DCT matrix generation How can I generate an NxM DCT matrix? In Matlab dctmtx function generates only square matrix (NxN).

type dctmtx

you will get:

function c = dctmtx(n)
[cc,rr] = meshgrid(0:n-1);
c = sqrt(2 / n) * cos(pi * (2*cc + 1) .* rr / (2 * n));
c(1,:) = c(1,:) / sqrt(2);

So it's not that hard to generate an non-square DCT matrix, just with minor revision:

function c = dctmtx1(n,m)
[cc,rr] = meshgrid(0:n-1, 0:m-1);
c = sqrt(2 / n) * cos(pi * (2*cc + 1) .* rr / (2 * n));
c(1,:) = c(1,:) / sqrt(2);


What does the convolution mean, what is the convolution philosophy I would like to know why the convolution is necessary.

that is, who said that multiplying numbers with others and then adding them would tell us something?

If you could give me analogies without mathematics, it would be a great help for non-mathematicians.

I'm somewhat in the same situation as you @Jony. But I'm working currently on a project involving audio processing, and I had to read a lot about convolution and related topics in the last few days.

I'm not sure I gained an intuitive understanding of convolution--esp. convolution vs. correlation-- but, to quote another answer, "Convolution in the time domain is equivalent to multiplication in the frequency domain." So it becomes obvious it may have a usage in filtering applications.

Wikipedia hs also a couple of nice animation that shows you how the convolution of two functions is related to the area of their intersection:

Depending on your field of interest, this may give you some applications idea. As the last word, let me repeat I'm very new to this topic, so take this answer with a grain of salt--and I would gratefully stand to be corrected by more experienced users