Computing Autocorrelation via FFT I have seen several methods to calculate Autocorrelations using FFTs, and am confused about why they differ.

Zero-Pad it to double its original length.Take the FFT. Then replace all the coefficients with their square magnitude, except for the first, which you set to 0. Now take the IFFT. Divide every element by the first one. This gives you the autocorrelation. https://stackoverflow.com/questions/4583950/cepstral-analysis-for-pitch-detection

create a centered version of x by setting x_cent = x / mean(x);
pad x_cent at the end with entries of value 0 to get a new vector of length L = 2^ceil(log2(N));
run x_pad through a forward discrete fast fourier transform to get an L-vector z of complex values;
replace the entries in z with their norms (the norm of a complex number is the real number resulting of summing the squared real component and squared imaginary component).
run z through the inverse discrete FFT to produce an L-vector acov of (unadjusted) autocovariances;
trim acov to size N;
create a L-vector named mask consisting of N entries with value 1 followed by L-N entries with value 0;
compute the forward FFT of mask and put the result in the L-vector adj
to get adjusted autocovariance estimates, divide each entry acov[n] by norm(adj[n]), where norm is the complex norm defined above; and
to get autocorrelations, set acorr[n] = acov[n] / acov[0] (acov[0], the autocovariance at lag 0, is just the variance). https://lingpipe-blog.com/2012/06/08/autocorrelation-fft-kiss-eigen/

And of course there is simply the first method but without setting the first element in the frequency domain to zero (i.e. what you'd expect from the Wiener–Khinchin theorem) https://stackoverflow.com/questions/3949324/calculate-autocorrelation-using-fft-in-matlab.
Why do the first two methods differ from the naive implementation of Wiener–Khinchin theorem. Are there any references that may help me figure out which method is best (and maybe on fast convolutions in general)?

The first method zeros out the "DC" in the signal see. If you compute the autocorrelation of $sin(\pi x / k) + 1$ using matlab's xcorr, you do not get a negative correlation for the kth lag. This doesn't make sense as a statistical correlation so you need to remove/zero out the "DC" in the spectrum corresponding to the $+ 1$.
The second method does not double the length and so will alias the signal. Maybe it is correcting for the aliasing. But also there are different definitions of autocorrelation see and may also be calculating the statistical autocorrelation assuming stationarity.


Is it worth recording 4K resolution from a 1080p camera input? I work with this as a recording setup: Canon 5D3 into an ATEM MINI recorder into Ecamm Live (primarily for the sake of bypassing the 30 minute limit at present, not for live-streaming), but eventually uploaded to YouTube.
I've been recording to the 1080p setting, but I'm concerned that my files are too compressed to start with, giving fairly poor quality on YouTube after export.

Ecamm Live has an option to record to a 2160p 4K resolution, which of course gives larger files sizes with higher bitrates.
Is this going to improve the quality of the image being recorded, where the clean HDMI input from the 5D isn't changing, i.e. is the higher bitrate going to produce a more uncompressed image than just recording to 1080? Or is it simply making a pixel of the same quality take up the space of 4?
Perhaps would uploading the final version in 4K to YouTube help improve the image?
Or is YouTube's compression such that there's not much I can do?

I do not know how the choice between 1080p and 4K affects the recording quality in Ecamm live.
But it makes sense to convert a high-bitrate 1080p video to high-bitrate 4K video before uploading to YouTube. In my experience, the resolution increase will significantly improve YouTube video quality. This is because YouTube will use a significantly higher re-encoding bitrate at 4K compared to 1080p. Only a part of the re-encoding bitrate increase is consumed by an increase in the baseline bitrate demand when going from 1080p to 4K. Much of it will remain available for improving the video quality.
Here is as an indirect example a busy frame of a high-bitrate 1080p video converted to high-bitrate 4K video and uploaded to YouTube, and played back using 1080p and 4K settings:

The blocky artifacts in the 1080p stream are much more prominent.


With the same channel S/N, the normalized mean square error differs I am working on a research work on transmitting ECG signals over wireless body area network. The signal is affected by noise and small scale fading. The normalized mean square error(NMSE) is used to estimate the quality of the reconstructed signal at the receiver. When I run the same m file in matlab several times, I obtain different values for the NMSE, although the same channel signal to noise ratio is used.

I think this is reasonable, do you agree?
If so, how can I obtain a single value to represent the quality of the reconstructed signal at the receiver?
I want to plot a curve with the S/N at the $x$-axis and NMSE at the $y$-axis, how can I plot this curve if the value of NMSE differ with different runs of the same m file with the same value of S/N?

The problem you describe is likely due to the randomness of your noise: Even though the S/N ratio remains the same for each run, the actual realization of the noise is different in each run. Hence, the NMSE is different for each run.

The standard procedure here is to run the algorithm several times for a single S/N ratio, measuring all the obtained NMSEs. Then, you calculate the mean of the measurements to find the overall NMSE for a given S/N ratio. Additionally, you can draw errorbars (MATLAB errorbar) around each points, which indicate the standard deviation of each measurement.


inter symbol interference vs digital modulation I am a little confused  after my lecture today. My professor started off with explaining what inter symbol interference (ISI) is and he ended the lecture talking about the different types of digital modulations out there.

I am now confused with how these two actually come together, or is it two different topics that have no connection?

Is ISI about sending digital signals in a wire? Does ISI come into play when doing lets say QPSK? I thought the symbol signal sent in QPSK is just a sin signal. or it uses this ISI?

ISI is an (typically undersired) effect of the channel – that can be a cable, an optical fiber, free space or the air you use to transport your wave.

The different kinds of constellations are a method of putting data in a waveform.

So, they are connected: The modulations define which symbols your system uses. ISI is what happens to the symbols when you rapidly send symbols one after the other.

Hence, your receiver needs to deal with ISI, and it needs to understand how to get the data back from the symbols it received.


A query on the non-uniform quantization I have read that non-uniform quantization boosts the smaller amplitude signals by a large amount. However the larger amplitude signals receive a small gain. As shown in the below diagram (Compressor Input and Output).

The Input - Output Characteristic diagram of compressor is given as below.

My query is as seen from the Input - Characteristic diagram, the step size is very less for smaller amplitude signal and hence the gain has to be less. On the other hand the step size is more for the large amplitude signal. Hence the gain has to be more. But the compressor input and output are in reverse way. Can you please help where I was going wrong?

To make it more clear, I suppose your question is

  Why it is said that the compressor gain at low input amplitudes is higher, while the step size of a nonuniform quantizer is small in that region. Similarly, Why it is said that the gain of the compressor is higher for high input amplitudes, while the step size is larger for those inputs.

First, notice that nonuniform quantization is used to improve the average SNR by adapting the quantizer to the source probability density function. Of course, when the source pdf is not uniform (e.g. when it is Gaussian), then uniform quantization is not a good choice. To better understand it, assume that you want to quantize a Gaussian source into $K$ levels. Let's first assume we have a uniform quantizer. One way to reduce the distortion caused by qunantization is to increase the number of steps ($K$). Increasing $K$ means the step size is reduced which eventually leads to a higher SNR (finer quantization). Increasing the number of levels means that the total number of bits to represent each sample is increased. For instance if we have 256 levels we need 8 bits to represent each level. But for 1024 levels we will require 10 bits.

Now in the above example assume that we want to work with a fixed $K$ (e.g. 256) and improve the SNR. To do this, notice that a Gaussian source is more likely to produce samples around the mean. So for example, a zero-mean Gaussian source is very likely to generate values close to zero. So we allocate smaller step size in the area around zero, and instead make the step size larger at the tail of the pdf (where it is less likely to happen). So technically, we quantize with a smaller step size on average (but without increasing the number of levels).

Any scalar nonuniform quantizer can be implemented by using a chain of

1- A compressor: a memoryless monotonic nonlinearity $G$

2- A uniform quantizer

3- An expander: inverse of nonlinearity $G$

In your example, the nonlinearity has a large gain at the origin (the characteristic function is steep) while it causes some sort of saturation effect at large amplitudes. But this gain is applied on the input signal not on the quantizer. Still the result should be quantized by a uniform quantizer afterwards. Although it is a uniform quantizer, but the preliminary change caused by the nonlinearity $G$ and then $G^{-1}$ makes it seem as if different step sizes are used: The low amplitude bins are magnified first, then there is a uniform grid, then they are compressed. They seem as they are quantized by a small step size. Something similar also applies to the high amplitude bins, except they are compressed first, then they are expanded. So the overall result is equivalent to the effect of a large step size.


Gaussian Blur In MATLAB and Connection to Image Resolution I am calling fspecial() in MATLAB to create a Guassian blur of an image.

For example,

%Create Gaussian Filter
G = fspecial('gaussian', [5 5], 2);

% Blur Image
blurredImage = imfilter(nonBlurredImage,G,'same')

The parameter [5 5] is hsize which defines the extent to which the Gaussian filter is applied ... which according to user "Try Hard" ... "In the case of a gaussian filter, the intensity at each pixel around the central one is weighted according to a gaussian function prior to performing a box average over the region."

The parameter 2 is sigma ... with units of pixels.

Say the size of nonBlurredImage is 200x200. And that the side of each pixel represents a length of 0.25 mm. Then the image could resolve objects with a "resolution" of 1 line pair/mm without aliasing.

If I pass G = fspecial('gaussian', [5 5], 1) to the filter function, and then produce the blurred image. What would be the "resolution" of this image? How does it relate to hsize? If I instead used [10 10] or [1 1], how would that affect things?

Thanks!

When applying Gaussian Blur on an image you should care for the STD only.
The rest should be set as following:

STD_TO_RADIUS_FACTOR = 5;
kernelRadius = ceil(STD_TO_RADIUS_FACTOR * kernelStd);
kernelLength = (2 * kernelRadius ) + 1;
mGaussianKernel = fspecial('gaussian', [kernelLength, kernelLength], kernelRadius);

Now, regarding the resolution of the image, this is tricky since it has many parameters in it.
If you want to measure it in Cycles in Spatial (As DFT), you can see the spectrum properties (Cut Off Frequency) of a Gaussian Kernel.

Update

The STD_TO_RADIUS_FACTOR should be determined by the balance between accuracy and performance.
The above is basically a truncated version of the Gaussian Kernel (Which is infinite). The balance sets how good is the approximation vs. how long the Convolution Kernel is (The longer the higher computational resources it requires).


What Is the Transfer Function of a Moving Average (FIR Filter)? To make post-processing easier, I export scope measurements as CSV files, which are then post-processed (mostly in Microsoft Excel, which is not the best tool for the job, but it is all I have at my disposal).

One of the post processing used is a simple moving average, and I was wondering how it does transform the source measurement. I expect it to behave similarly to a low-pass filter, but I guess it is not the same. Hence, my question.

I am specifically looking for a Bode plot. Thanks in advance.

The frequency response of the moving average is called the asinc or psinc, the aliased sinc or periodic sinc (sinc for cardinal sine), or the Dirichlet function.

Since the sum of the moving average filter coefficients is equal to one, it preserves constant signal, hence is somewhat lowpass. When the length $N=2$, the shape in the frequency domain is shaped like a cosine, hence decreasing.

Otherwise, it ripplles in the frequency domain. When the length is even, the amplification at Nyquist is $0$, when odd, of $1/N$, as shown below:

So it's globally lowpass, not the best, but one of the fastest to compute, allowing fast recursive computations.


About sampling and wireless networks Wireless networks operate on a high frequency, say some several gigahertz. Are the signals really sampled at that high frequencies? If I got it correctly there's some modulating/demodulating in signal transmit, does that effect the situation somehow?

What you describe is a direct-sampling architecture. Those are practically available for radio amateurs using frequencies up to 30MHz. They are also actively researched to open up their usage for much higher frequencies. Many interesting topics (e.g. using the aliasing effect to your advantage and undersample the signal) are investigated here and the possibilities of these receivers are large. However, it will probably take a long time until those are practically available for the WiFi frequency range.

However, your WiFi transceiver is most likely a direct-conversion architecture. In the digital domain of it, the data-carrying signal is generated. In case of 802.11n it is an OFDM modulated signal with up to 20/40MHz bandwidth. This signal is centered around 0Hz and hence cannot be transmitted as it is. This signal is called the baseband signal.

After Digital/Analog-conversion, in a stage called a mixer, this signal is then upconverted to the final carrier frequency around 2.4GHz. In the spectrum, the signal is literally shifted from a center frequency of 0Hz to a center frequency of 2.4GHz. At this point, the signal is called a passband (or bandpass) signal. After amplification, this signal is then transmitted via the antenna.

On the receiver side, the same happens the other way around. There is another mixer that downconverts the passband signal from 2.4GHz center frequency to a baseband signal at 0Hz center frequency. At this point, the signal has just a bandwidth of 20/40MHz and can easily be sampled in an A/D converter.

What follows is demodulation, error correction and other digital processing to eventually show you pictures of cute kittens in your browser.


Why do we integrate the modulating function in FM? I've been reading about FM and I keep coming to an expression, where the instantaneous phase

$$\Theta=2\pi f_it+\phi\text, \tag1$$ where
$$f_i=f_c+k_fm(t)\text,\tag2$$
($f_c$ being the base carrier frequency and $k_f$ being some proportionality constant) is written as some sort of integral expression, like here:

$$\Theta_i\{t\}=2\pi f_ct+2\pi k_f\int_0^t m\{t\}\,dt$$

No clear explanation is given, apart from:

  The instantaneous phase of the modulated wave at any instant can be obtained by substituting 12.6 into 12.3 and integrating to get [the integral expression]

(12.3 and 12.6 being my equations (1) and (2) respectively).

Why do we need to integrate the modulator signal?

Why is the following not an accurate expression for the phase: $\Theta =2\pi f_ct+2\pi k_fm(t)t+\phi$, (simply plugging equation (2) into equation (1))?

You need to integrate the modulating signal because frequency is the time derivative of phase. Therefore, the typical relationship from introductory calculus holds:

$$
\phi_i(t) = \int_{-\infty}^{t} \frac{d\phi_i(\tau)}{d\tau} d\tau = \int_{-\infty}^{t} 2\pi f_i(\tau) d\tau
$$

For causal signals, the lower limit on the integral changes to zero:

$$
\phi_i(t) = \int_{0}^{t} 2\pi f_i(\tau) d\tau
$$

$$
\phi_i(t) = \int_{0}^{t} 2\pi\left(f_c + k_fm(\tau)\right) d\tau
$$

$$
\phi_i(t) = 2\pi f_ct +\int_{0}^{t} 2\pi k_f m(\tau) d\tau
$$

which is the answer you were trying to get to.


Why does pulse-width modulation work? When we apply pulse-width modulation (PWM) on controlling the brightness of a LED or the speed of a DC motor, is it the lag in our eyes that is responsible for the feeling of a continuous powering or the inertia of the devices?

Depends on what you are controlling. For DC-motors it is the inertia of the device that acts as a low-pass filter of the PWM modulated signal resulting in a continuous motion. For most LEDs it is the human eyes that do the apparent low-pass filtering.

If the PWM-frequency is not very high you can actually see this by moving your head from left to right while looking at the PWM modulated LED. The light will appear on and off while the LED is moving through your field of vision (this can get very annoying with LED car taillights).


Expectation and Auto-correlation of an Output from a Matched Filter. [see the attached image]

I recognize the filter in the problem is a matched filter. But I don't understand how this filter actually work, and I am not quite sure about how to calculate the expectations and the auto-correlation in this case. Any help is greatly appropriated.

First you should find $Y(t)$. In this case, it can be calculated as:

$$Y(t)=X(t)*h(t)=\int_{-\infty}^{\infty}h(t)X(t-\tau)\ \mathrm{d}\tau$$

But you have more information about $X(t)$ and $h(t)$:

$$Y(t)=\int_{-\infty}^{\infty}S(-t)[S(t-\tau)+N(t-\tau)]\ \mathrm{d}\tau$$

If we assume that the integral converges (which makes sense, because if not $Y(t)$ would not exist), then we can do the following:

$$\begin{align}
\mathbb{E}[Y(t)]&=
\mathbb{E}\left[\int_{-\infty}^{\infty}S(-t)[S(t-\tau)+N(t-\tau)]\ \mathrm{d}\tau\right]\\
&=\int_{-\infty}^{\infty}\mathbb{E}\left[S(-t)S(t-\tau)+S(-t)N(t-\tau)\right]\ \mathrm{d}\tau
\end{align}
$$

But $S(t)$ is deterministic, so:

$$\mathbb{E}[Y(t)]= \int_{-\infty}^{\infty}S(-t)S(t-\tau)+S(-t)\mathbb{E}[N(t-\tau)]\ \mathrm{d}\tau$$

The problem states that the mean of the noise process is zero, thus:

$$\mathbb{E}[Y(t)]= \int_{-\infty}^{\infty}S(-t)S(t-\tau)\ \mathrm{d}\tau$$

So, the expectation of $Y(t)$ is the convolution between $S(t)$ and its reversed version (also known as autocorrelation of $S(t)$):

$$\mathbb{E}[Y(t)]=S(t)*S(-t)$$

Because $S(t)$ is known, then $\mathbb{E}[Y(t)]$ is, too.

To find the autocorrelation, you have to do something similar to this procedure. I'll leave that to you.


Opinion about this aproach on Machine Learning I have a question for you that maybe you could give me a clue.

I'm trying to make a CBIR (Content-Based Image Retrieval), so I query an image and I get the most similar, not classified within a database of around 2000 images. I'm trying to extract features and matching them by euclidean distance but I was going to try to train them with the framework called Caffe, for deep learning (neural networks).

So, what I've done is taking 10 samples per every image by 3d projection (so I take the image 635.jpg and I make 635_1.jpg and so on to 635_10.jpg) and I use the original image 635.jpg as validation data.
So if there are 2000 images, I use 20.000 images for training and 2.000 images for validation into 2.000 clases... do you think it's a good idea?

Thank you in advance.

Although I've been trying, I wasn't still able to train on Caffe due to the high quantity of errors in the dependencies on Mac OSx. I will switch to Ubuntu and see what's wrong.

This is a fine idea and pretty standard. However, your problem will be the augmentation of the data. If you render the image from many viewpoints/rotations etc. and just use a plain background, the network would include that information into the training stage e.g. Like a shape cue. So, you have to use different backgrounds. This is of course scene specific. However, you might be able to narrow down your scenes so that you could randomly generate/perturb different background structures.

Also note that, you might be able to scale up to thousands of classes, but getting to millions could become a headache, since training this network would be very tough.

You might as well go with standard indexing / hashing approaches, which already perform pretty well.


Fill matlab "ellip" function through transfer function magnitude response The matlab "ellip" function can be used to design the unquantised coefficient set.

  From matlab website:

  $ [b,a]=ellip(n,Rp,Rs,Wp) , $

  n: order of filter

  Wp: normalized passband edge frequency

  Rp: ripple in passband (db)

  Rs: db down from he peak value in the passband

If i have an 8th order IIR filter with the following transfer function magnitude response:

\begin{align*}
 0\,\mathrm{db} \pm 1.5 \,\mathrm{db} & \quad \text{for } \quad 0.2<|v|<0.3 \cr
 \lt −60\,\mathrm{db} & \quad\text{for } \quad |v|<0.14\text{ and } |v|>0.36
\end{align*}

how i can fill the matlab function with the suitable values?

The quickest way would be

[b,a] = ellip(8,1.5,60,[.2 .3]);

This designs a bandpass from .2 to .3 with 1.5 dB ripple and 60 db stop band attenuation. You don't get to pick the edges of the stop band since you already pre-selected the order. You can do one or the other but not both at the same time. However the 8th order filter is more than enough to meet your stop band requirements. In fact the low band stops at 0.195 and the high stop band starts at 0.307.

Another wrinkle is that the ripple actually only goes down, i.e. passband is not between -1.5dB and +1.5 dB but only between -1.5dB and 0 dB. In order to work aound this (if you care) you can double the ripple to 3 dB and than simply add 1.5 dB gain. You would also have to increase the stop band attenuation to 61.5dB to compensate for the overall gain change. The funcion $ellipord$ allows to calculate the required order to achieve a certain stop band, so the correctest way would be

%% Ellip filter design example
[n,wp] = ellipord([.2 .3],[.14 .36],3,61.5);
[b,a] = ellip(n,3,61.5,wp);
% add 1.5 dB of gain
b = b*10.^(1.5/20);
[H,F]= freqz(b,a,1000,2);
plot(F,10*log10(H.*conj(H)));
set(gca,'ylim',[-70,3]);
grid('on');

So it turns out a 5th order elliptic filter can do the job just fine.


What is the no overload region of a quantizer? I was watching a video on quantization noise and the instructor kept referring to the "no overload" region of a quantizer? Could anyone please explain what that term means?

It is stated in chapter 5 (Scalar Quantization I: Structure and Performance) of the book "Vector Quantization and Signal Compression" by Gersho-Gray that

  Associated with every $N$ point quantizer is a partition of the real line
  $\mathcal R$ into $N$ cells or atoms $R_i$, for $i = 1,2,\cdots, N$. The $i$th cell is given by $\mathcal R_i=\{x \in \mathcal R: Q(x) = y_i\} \equiv Q^{-1}(y_i)$, the inverse image of $y_i$ under $Q$. It follows from this definition that $\bigcup_i \mathcal R_i = \mathcal R$ and $R_i\bigcap R_j = \emptyset$ for $i \neq j$. A cell that is unbounded is called an overload cell. Each bounded cell is called a granular cell. Together all of the overload (granular) cells are called the overload region (granular region).

So in other words, the union of bounded regions (i.e. granular region) where saturation does not occur is referred to as the "no overload" region and is identified in the following picture


Match (OCR?) a digital image of a product label to similar Excuse my ignorance...

Given a digital image of a label I'd like to

Requirement: Match it to the closest pre-existing image in a datastore
Optional: OCR it, if possible (only necessary if step 1 fails)

I'm a programmer but not an engineer. I have no idea where to begin or where to look into solving this problem. I don't even know the right terminology; "Image Processing" is best I could come up with.

Starting here, on DSP (beta); apologies in advance to the forum police if this type of inquiry is not considered to be 'in scope'. Before you vote me off the 'island' hopefully someone can comment with a better site/forum for my education.

First, you will need an image processing algorithm to detect local features of the digital images you have. 1 and 2 provide good reviews of such algorithms.:

After extracting the local features, a matching process between the extracted features and your set of labels is used to link the labels with acquired images. The second references goes over few matching techniques and performance metrics.


MPEG Layer I: use of scalefactor and quantizers I'm studying MPEG compression on the book "Introduction to Data Compression" by Khalid Sayood.
In the figure there is a paragraph piece about the subject.

I didn't understand a few things:

use and how the scalefactor works
use and how the 14 different quantizers work.

That are the parts highlighted in orange.

The scalefactor is a value that is used to decrease (or increase) the range of values, right?

I try to repeat what I understand:

Input signal is PCM samples over time
These samples are transformed into the frequency domain (I suppose with an FFT)
Samples are subdivided into frames of 384 samples each
Each frame is processed individually and subdivided into 32 equal frequency bands using a filter bank. Each subfield therefore contains 12 samples
encoder processes each sub-band individually and through a psychoacoustic analysis determine the number of bits to allocate for each quantized sample and possible quantization levels (linear quantization)
encode the samples and form the package
etc.

I know that point 5 is very superficial because I didn't understand the parts related to the scalefactor and the 14 different quantizers.
Can anyone help me?

Somewhere I've read the quantization steps are the same for each sub-band and they are 2 dB each (linear quantization), so what's the scalefactor meaning?
I think the scalefactor is a kind of exponent in the case of non linear and floating point quantization.
So MPEG layer I use a linear or not lineare quantization?

That's what I understand. Why does the book mention 14 quantizers?

I read (quickly) the book you've been advising me (Introduction to Digital Audio Coding and Standards from Bosi).
This is a piece that treats the part I didn't understand.

For each sub-band, the scalefactor (stored on 6 bits) is calculated from the 63 available and a quantizer (on 4 bits) of the 14 available.

The scalefactor is the largest sample of a sub-band.
Once I've calculated it I use to divide each sample, obtaining 12 samples whose value is in the range [0, 1].

The quantizer has a numeric value that expresses the number of bits I use to encode each sample.

This is my usual scheming:

Why do I need to normalized the samples?

Let's get back to the official MPEG standard document and avoid your confusing literature for a moment!

So, you get audio as time samples. You use a very specifically crafted 32-band filterbank, which takes these input sample stream and generates 32 sample streams of each 1/32 of the original rate (that's the "time/frequency mapping" in that diagram you show, and the QMF in this diagram). All these streams are time-domain!¹

Now, you take a frame – that is, 12 "parallel" samples from each of these 32 streams. The codec's job is to figure out how "important" they are, and only use as much bits as needed to represent these.

The "importance" step is done by the psychoacoustic model (which might do arbitrary complex things, but that's not the point of discussion here), which will give us some info of the type "use $n_0$ bits for the first subband, $n_1$ for the second,…".

The Quantizer's job is now to take the 12 bits of "its" subband and save them with the $n$ bits depth it got. That's pretty straightforward. You just take the highest value in your 12 samples, scale it so that it is $2^n$, scale the other 11 samples as well, round them to integers, and save the scalefactor for later reconstruction.

You do that to all 32 subbands, and send the resulting streams of bits together with the scalefactors as the data part of an MPEG 1 Layer I frame.

At the decoder, you multiply with 1/scalefactor, then you transform the 32 subbands back into a single time signal.

Your confusion stems from the fact that you took the intro to that chapter as a kind of "step before" encoding happens, whereas it just was a slightly misleading (and claiming the samples were transformed to frequency domain) overview.

Also, the author of your book is wrong about the usage of the FFT in MPEG 1.

Be very precise about what you're exactly considering. We're here talking about MPEG 1 Layer I, not Layer III. In layer III ("MP3"), there's a kind-of-a-frequency-domain-transform happening on the subband samples (a modified Discrete Cosine Transform). That might have contributed to confusion.

¹ the "frequency mapping" is actually the mapping of signals to the frequency-disjunkt subbands. There's a bit of a mathematical caveat to that – the QMF filterbank is pretty close to being a wavelet transform, but I really wouldn't call it a time-frequency transformation. The original standard wouldn't want to do that, either, so they intentionally used the word "mapping". By the way, you'll find more info in the standard, which you can find by following the picture source that the author of your list (presumably) has in its table of figures.


What are the units of a .wav file? I am working on calibrating my recording system. I have a signal generator playing a 500mVpp sine wave signal into my utrack24 recorder. I am trying to determine what the units of the recorded utrack .wav file is so I can compare that to the input signal. I have read the .wav file into matlab with audioread. I can use either the normalized version of the file or the 'native' version; I'm not sure which is better.

Once an audio signal has been digitalized, it has no units and it is treated as a number.

You could use dBFS (decibels relative to full scale) as unit, though. The value of 0 dBFS represents the maximum digital signal level.


Library/strategy for rotate & scale image alignment/registration I have two raster images of maps, A and B.  I know that B contains A but I do not know the scale or rotation.  There should not be any other distortions between the image but there may be noise.  Please could anyone recommend an approach or software library which would provide me with the scale and rotation transformations to align A onto B, thus allowing me to transform arbitrary points between the coordinate spaces?

Left: A: source, B: destination.

Right: Location of A on B

Given the high contrast, plane sharing and 'edgyness' of the images I would hope it would be fairly trivial compared to most image analysis tasks.  There are similar questions on here but they seem to be doing something more complex.

One option is keep all the edge points, or to identify key points (junctions, corner) on both images, and then perform a (semi-)rigid matching, such as Coherent Point Drift, which is robust to missing points.

Further, you could split the biggest images into smaller overlapping patches, and adapt the template to the different patches.


Why need an low pass filter after up-sampling? I understand why there need an low pass filter before down-sampling, the sample frequency need to be at least twice of the max frequency signal, else there will be alias issue.

But base on the Wiki page, we need an lowpass filter after up-sampling.

  Upsampling requires a lowpass filter after increasing the data rate,
  and downsampling requires a lowpass filter before decimation.

I don't understand that, can anyone explain it?

Let's say you've sampled an analog signal $x(t)$ with spectrum $X(\omega)$ at rate $1/T$ that is high enough to satisfy the sampling theorem. The spectrum (i.e. the discrete time Fourier transform (DTFT)) of the sampled signal $x_{1,k}$ will be a periodic repetition of $X(\omega)$. The repetition period is $1/T$.

Now you sample the same signal with a higher rate $L/T,\, L\in\mathbb N$ yielding $x_{L,k}$. Again the spectrum will be a periodic repetition of $X(\omega)$ but this time the repetition period is $L/T$, so the spectral images have greater frequency distance than before.

The task of upsampling consists in calculating $x_{L,k}$ from $x_k$. First, $L-1$ zeros are inserted after every sample of $x_k$. Actually this just changes the basic frequency support of the DTFT to $-L/(2T)\ldots L/(2T)$ containing $L$ copies of the original (analog) spectrum. Therefore the unwanted copies are filtered out with a lowpass filter so that only the original spectrum in range $-1/(2T)\ldots 1/(2T)$ remains. This is identical to $x_{L,K}$

The above steps are quite well explained in the figure of the Wiki article you quoted (in the same order).


Pyramid vs Scale-space in Focus Stacking I'm working on a focus stacking program, I wrote one a few years ago which used a very simple method, basically taking a DoG approach, but only operating on the single full resolution image. The resulting images had a soft focus look, kind of dreamy.

This time I researched further and found papers describing image pyramids, essentially doing the same thing as before except on multiple spatial frequencies, then recombining them with the maximum components which I imagine will solve the soft focus look.

But learning about pyramid based approaches led me to scale-space representation and multi-resolution analysis, and I can't tell if I should pursue one of those methods for possibly better results? (And MRA led me to wavelet transforms which I also don't see how I would implement.)

Furthermore, I don't really understand how either would even be implemented for focus stacking, but my impression of scale-space is that in practice I would just end up building a pyramid anyway, but that doesn't seem right. So am I missing something about that?

And does it seem like there is anything to be gained by using a newer or more sophisticated approach than image pyramids for focus stacking?

(I used Python last time but am using Apple's CoreImage framework this time, so I might be a bit more limited now, maybe MRA isn't even an option.)

Thanks in advance!

Your question relates to techniques with different names. The  first one is related to micrographs taken at different (aligned) focal positions combined  into a focus image. This Extended Depth of Field is based on complex wavelet-based methods.

To answer your second question on difference between MRA, pyramids, etc.: most pyramids and wavelet-based methods for 2D image processing are instance of multiscale/multirate//multiorientations  decompositions based on certain scale-space discretizations (a review paper on these issues is A Panorama on Multiscale Geometric Representations, Signal Processing, 2011). With scale-space in general, you do not need to downsample, and end up with a stack, not a pyramid.

Since one selects maximum coefficients across scale, band, or orientation, my opinion is that a complex pyramid, a complex (or dual-tree) wavelet transformation, with a little redundancy and some shift or rotation invariance could be useful, at the cost of more painful implementation. Standard real MRA wavelets are limited in orientations, full scale-space might be quite involved.

Maximum coefficient selection is probably  a poor's man saliency detection, which is another track to follow based on the specificity of your images.

More general terms for your search are multifocus image fusion, multiscale or multi-resolution images fusion, even multichannel. A few references:

Performance comparison of different multi-resolution transforms for image fusion, Information Fusion, 2011
Multifocus image fusion and restoration with sparse representation, IEEE Transactions onInstrumentation and Measurement, 2010
Multifocus image fusion using the log-Gabor transform and a Multisize Windows technique, Information Fusion, 2009

You can as well find inspiration in multiscale versions of High-dynamic-range (HDR) imaging, for instance:

Local Edge-Preserving Multiscale Decomposition for High Dynamic Range Image Tone Mapping, IEEE Transactions on Image Processing, 2012
Multiscale Processing with HDRWaveletTransform


Is DFE equalizer always blind or can we use reference signal? I have made a blind DFE Equalization, I mean decision directed .
But I am not sure how to make an implementation based on reference signal known by receiver and emitter. In the back forward filter adaptation how do I take this into consideration?

My second equation is to know if the updating rule for the RLS is the rigth one ?

Here is my rls tracking DFE matlab code

`

    clearvars;
    close all; clc;
    MQAM = 2;
    N = 1e2;
    L = 3;
    LFF = 8;
    LFB = 8;
    target = 1e-6;
    error_min = 30;
    trans_min = 1e6;
    trans_max = error_min/target;
    BER=0;
    SNR=0;
    iSNR=1;
    endSNR = 16;
    stepSNR = 2;
    h = [1 0.5*exp(1i*pi/6)  0.1*exp(-1i*pi/8)]';
    eq = comm.DecisionFeedbackEqualizer('Algorithm', 'RLS', ...
               "AdaptAfterTraining",false, 'NumFeedbackTaps', LFB, ...
               'NumForwardTaps', LFF); eq.ReferenceTap = 1;
    BER_mtlb = 0;
    handle = waitbar(0, 'Please wait...');

    p=LFF; lamda=0.98; sigma=1; P1=(sigma^-1)*eye(LFF);

    P2=(sigma^-1)*eye(LFB);

    while SNR(iSNR)<= endSNR
         error_num=0;trans_num=0;  error_mtlb = 0;
         while  (error_num <= error_min) && (trans_num <= trans_min)
             MemFF = zeros(LFF,1); MemBW = MemFF; hFF = zeros(LFF,1);
             hFB=hFF; hFF(ceil((LFF-1)/2)) = 1;

             s = randi([0 1],N,1);
             smod = qammod(s,MQAM,'gray','InputType','bit');
             xin = awgn(filter(h, 1, smod), SNR(iSNR), 'measured');

             out = zeros(N,1); dec = zeros(size(s)); sdec = zeros(size(s));
             smodCplx = complex(smod);
             [y, err, weigths] = eq(xin, smodCplx);
             rxx = qamdemod(y, MQAM, 'gray', 'OutputType', 'bit');
             error_mtlb = error_mtlb + sum( rxx ~= s);
             %- RLS DFE Equalization loop block
             for n = 1: length(xin)
                 MemFF = [xin(n); MemFF(1:LFF-1)];
                 out(n) = transpose(hFF)*MemFF + transpose(hFB)*MemBW;
                 sdec(n) = qamdemod(out(n),MQAM,'gray','OutputType','bit');
                 dec(n) = qammod(sdec(n),MQAM,'gray','InputType','bit');

                 alpha(n) = dec(n) - out(n);% error of RLS

                 g1(:,n)=P1*MemFF*((lamda + MemFF'*P1*MemFF).^-1); % Gain
                 P1=(lamda^-1)*P1 - g1(:,n)*MemFF'*(lamda^-1)*P1;

                 g2(:,n)=P2*MemBW*((lamda + MemBW'*P2*MemBW).^-1); % Gain
                 P2=(lamda^-1)*P2 - g2(:,n)*MemBW'*(lamda^-1)*P2;

                 hFF = hFF  + alpha(n)*g1(:,n);
                 hFB = hFB + alpha(n)*g2(:,n);
                 MemBW = [dec(n); MemBW(1:LFB-1)];

             end

             M = min(length(s), length(sdec));
             Ryx = xcorr(1-2*sdec(1:M),1-2*s(1:M),'unbiased');
             [valMax,Delay] = max(abs(Ryx(M:M+99)));
             Delay = Delay-1;

             if (Ryx(M+Delay) < 0)
                 sdec = qamdemod(-out,MQAM,'gray','OutputType','bit');
             end
             sdec_synchronized = sdec(1+Delay:end);
             s_synchronized = s(1:length(sdec_synchronized));

             error_num = sum(s_synchronized~= sdec_synchronized) + error_num;
             trans_num = trans_num  + length(sdec_synchronized)  ;
         end % end this number of transmitted symbols

         BER(iSNR) = error_num/trans_num;
         BER_mtlb(iSNR) = error_mtlb/trans_num;
         waitbar(log10(BER(iSNR))/log10(target),handle);
         iSNR = iSNR+1;
        SNR(iSNR) = SNR(iSNR-1) + stepSNR;
 end
    SNR=SNR (1:length(BER)); berTheory =
    berawgn(SNR/log2(MQAM),'psk',2,'nondiff');  SNR =
    SNR(1:min(length(BER), length(BER_mtlb))); h = figure;
    semilogy(SNR,BER, 'b--'); hold on; semilogy(SNR, BER_mtlb, 'g--'); hold on;
    semilogy(SNR, berTheory, 'kx-'); legend('Simulation DFE', 'Theory');
    xh=xlabel('E_b/N_0 (dB)'); %yh = ylabel('Probability of bit error');
    xlabel('SNR (dB)'); ylabel('BER'); grid on;
    XMIN= 0; XMAX = 30; YMIN= 1e-7; YMAX= 1; axis([XMIN XMAX YMIN YMAX])

I can answer for adaptive filter containing feedback path alone as I have implemented it recently (both Blind DFE and with Reference Symbols), but it should work for DFE containing feed-forward path as well. In Blind DFE while computing the error, $e[n] = x[n] - \hat{x}[n]$, where $x[n]$ is the equalizer output, $\hat{x}[n]$ is the corresponding Hard Decision. But if you have reference symbol $r[n]$, you need not compare with hard decision output and $e[n] = x[n] - r[n]$.

If the reference symbols are placed in the beginning of data, then the filter adapts more reliably (delay and tap values) in the beginning itself, you can then use your Blind DFE for the remaining data part.

The example I have used comes from 802.11ad Single Carrier data packet where for each frame containing 512 symbols, the first 64 are known reference symbols and the rest 448 are random data symbols (BPSK/QPSK/16QAM). So for the first 64 symbols, I skipped the hard decision slicer and directly compared adaptive filter output with reference symbols. For the remaining 448 symbols, I used Blind DFE.


Independent Subspace Analysis and Topographic ICA Can anyone recommend the best references for independent subspace analysis (ISA). I have read 'Multidimensional Independent Component Analysis' by Cardoso and would like to know of any modern, solid (i.e. generally accepted to work) approaches.

If you have implemented ISA before how did it perform?

Read this paper about ISA: http://www.cs.helsinki.fi/u/ahyvarin/papers/ESANN06.pdf

and this one about topographic ICA: http://www.cs.helsinki.fi/u/ahyvarin/papers/NC01_TICA.pdf

Hyvärinen is the author of the FastICA algorithm.

PS.

By the way, you might want to look at ICASSO (http://research.ics.aalto.fi/ica/icasso/), if you are not already familiar with it.


Add echo as described by difference equation to audio signal (in MATLAB) A simple linear system is echo. It can be described by equation
$$
y[n] = x[n]+ k\,x[n−d],
$$ where $n$ represents sample index, $k$ represents an attenuation coefficient, and $d$ represents time lag.

How can I add echo to an input signal in MATLAB with convolution?

I'm new at MATLAB.

The impulse response $h[n]$ can be obtained by inspection:
$$
h[n] = \delta[n] + k \delta[n - d]
$$

In MATLAB notation this would be:

h = [1, zeros(1,d-1), k];

The output is computed by:

y = conv(x, h);

Note that the time sample indices for the vector y will be as follows, assuming that both x and h are indexed where the first element corresponds to time $n=0$:

ny = 0:(length(x) + length(h) - 2);


How to get the 10-5 most dominant frequency in an audio signal I have this wav file that I recorded using my laptop, It contains the sound of a G chord of a guitar. Since the guitar have 6 strings the sound should contain different frequencies I want to determine the top 10 most dominant frequency in the signal.

[X, Fs] = wavread('chordG.wav');
freq = fft(X);
plot(abs(freq)):

I don't know how to proceed after that it shows the graph of the frequency but I don't know how can I extract the values in the graph.

Your answer would be very much appreciated :)

There are research papers on the topic of automatic polyphonic chord estimation/transcription on MIREX (http://www.music-ir.org/mirex/wiki/2014:Audio_Chord_Estimation).  Also a book and a thesis by A. Klapuri.


two dimensional integration of a trigonometric function I am working on a detection problem which finally, I have to solve the following 2-D integral:
$$\int\limits_{a}^{b} \int\limits_{c}^{d}e^{A\sin(x)\cos(y-B)}\, \mathrm{d}x \, \mathrm{d}y \ ,$$
where $a,b,c,d$ ,$A$ and $B$ are constants. Does anyone can help me with the analytical solution for this integral?

Since, the Taylor series expansion for $e^z$ about $z=0$ is
$$e^z = \sum_{n=0}^{\infty}\dfrac{z^n}{n!}$$
then, ignoring the question of convergence, you can say
$$\begin{align*} \int_a^b \int_c^d e^{A\sin(x)\cos(y-B)}dx dy &= \int_a^b \int_c^d  \sum_{n=0}^{\infty}\dfrac{(A\cos(y-B))^n}{n!}\sin^n x \; dx \, dy\\
\\
&= \int_a^b   \sum_{n=0}^{\infty}\dfrac{(A\cos(y-B))^n}{n!}\int_c^d\sin^n x \; dx \, dy\\
\\
&=   \sum_{n=0}^{\infty}\dfrac{A^n}{n!}\int_a^b \cos^n(y-B)\int_c^d\sin^n x \; dx \, dy\\
\\
&=   \sum_{n=0}^{\infty}\dfrac{A^n}{n!}\left(\int_{a-B}^{b-B} \cos^n u\;du\right)\left(\int_c^d\sin^n x \; dx \right)\\
\\
&=  \left(u\big|_{a-B}^{b-B}\right)\left(x\big|_c^d\right) + A\left(\sin u\big|_{a-B}^{b-B}\right)\left(-\cos x\big|_c^d\right) \\
&\quad + \sum_{n=2}^{\infty}\dfrac{A^n}{n!}\left(\int_{a-B}^{b-B} \cos^n u\;du\right)\left(\int_c^d\sin^n x \; dx \right)\\
\end{align*}$$
Using integration by parts, the reduction formulas for the integrals of $\sin^nx$ and $\cos^nu$ work out to be
$$\int \sin^nx \; dx = -\dfrac{1}{n}\sin^{n-1} x \cos x +\dfrac{n-1}{n}\int \sin^{n-2} x \;dx$$
and
$$\int \cos^nu \;du = \dfrac{1}{n}\cos^{n-1} u \sin u +\dfrac{n-1}{n}\int \cos^{n-2} u\; du$$
For $|A|\gg 0$, you could use the Taylor series expansion of $e^z$ about $z=A$ for a series that converges more rapidly for that value of $A$.


What is the Concept of MATLAB Function Polynomial Interpolation? I am curious about MATLAB function.
Can you tell me why does not it do good approach to using of MATLAB function except for speed reason?

Is there any other reason to using polynomial utility?

It is basically an approach choice.
Inside the math is identical.

Usually, when doing Least Squares curve fitting, you're not looking for the Polynomial coefficients but a scaled version of them.
For instance, if you try to estimate Range, Velocity and Acceleration from Range Measurements using LS you need to scale the Polynomial Coefficients according to the dynamic model.

By the way, most of the times, you can do the LS estimation using $ {\left( {H}^{T} H \right)}^{-1} {H}^{T} $ which would be faster than pinv.
The only case it won't work is when $ \left( {H}^{T} H \right) $ isn't invertible (Rarely occurs in real world measurements).


Which one of these signals will be least distorted by this communication channel So I was solving a problem that resulted in a communication channel that heavily distorts high frequency but does not significantly affect low frequency and I got stuck at figuring which one of these signals would be the least distorted. I know I have to find the one with the lowest frequency but I'm not completely sure how. These are the signals:

$$x_1(t) = \delta(t);\quad x_2(t) = 5; \quad x_3(t) = 10e^{j1000t};\quad x_4(t) = \displaystyle\frac 1t$$

I just need to know how to find the frequency of one of the signals

Normally, a signal is composed of several frequencies, which strengths are given by the Fourier Transform of the signal.

$x_1(t)=\delta(t) \Leftrightarrow \mathcal{F}(x_1(t))=C$, i.e. the Fourier transform is a constant (value of constant depends on the actual normalization of the Fourier Transform). Hence, it has all frequencies contained, with equal power.
$x_2(t)=5 \Leftrightarrow \mathcal{F}(x_2(t))=\delta(f)$, i.e. the Fourier Transform of a constant is a Dirac impulse at frequency zero. Hence, this signal has only one frequency component, the zero-frequency.
$x_3(t)=10e^{j1000t} \Leftrightarrow \mathcal{F}(x_3(t))=C\delta(f-1000/2\pi)$ i.e. the Fourier Transform of a complex exponential is a Dirac impulse at the position of the frequency of the exponential. Hence, it has only one component, at frequency $1000/2\pi$
$x_4(t)=1/t \Leftrightarrow \mathcal{F}(x_4(t))=C\text{sign}(f)$, i.e. the Fourier transform of the inverse is a sign function. Hence, it contains all frequencies.

Accordingly, the answer to your question would be signal 2, as it only contains the DC frequency ($f=0$).


Directivity Pattern of microphone array Directivity pattern of a microphone array is fourier transfrom of  microphone weight values. By modifying the amplitude weights, we can modify the shape of the directivity pattern.

I can really see this when all microphone weights are equally valued.
Fourier transform of a rectangle is a sinc func like picture below:

Now I believe if I apply a sinc function for my microphone array's element values I should get a rectangular shape in my directivity pattern.

I have a simulation program of my own. Which works when microphone values are equally valued . But I am not getting what I expect(a rectangular pattern) when I give sinc values to microphone array elements.

Am I missing something teorically or it is just my program bugged?

Edit: As requested in comments, weights that I apply can be seen in the graph below. 40 is the microphone count of my array. Sorry for disgusting looking graph it is the output of my C++ program.

As I said in the comment, you sound like you are doing things correctly.

Below are two plots generated by the scilabcode at the end.

The first plot shows two different weightings applied to the sensors.  The second plot shows the associated beam patterns.  While the $\tt sinc$ one is not "uniform" it is a fair approximation.

EDIT

Changing the code below to cater for an even number of sensors:

NOn2 = 20;
N = 2*NOn2;
w_uniform = ones(1,N);
omega = 2*%pi*([0:N-1] - N/2)/N + 0.0001;
w_sinc = [sin(N/15*omega)./(omega) ]

Then I get the plots below.

The relationship between the weights and the beam pattern is like a Fourier transform:  the broader the sinc, the narrower the beam pattern. The narrower the sinc, the wider the beam pattern.

CODE ONLY BELOW

// http://www.astron.nl/other/workshop/MCCT/MondayPatel.pdf
lambda = 1;
d = lambda/3;
theta = [-%pi/2:0.01:%pi/2];

phase_adjacent = 2*%pi*d/lambda*sin(theta);

NOn2 = 7;
N = 2*NOn2 + 1;
w_uniform = ones(1,N);
omega = 2*%pi*[-NOn2:NOn2]/N + 0.0001;
w_sinc = [sin(N/4*omega)./(omega) ]

e_to_phi = exp(-%i*phase_adjacent.'*[0:N-1]);

r_uniform = e_to_phi*w_uniform';
r_sinc = e_to_phi*w_sinc';

figure(0);
clf;
plot(theta/%pi*180,abs(r_uniform));
plot(theta/%pi*180,abs(r_sinc),'r.');
title('Uniform weighted (blue solid) and sinc weighted (red dotted) array resposnes');

figure(1);
clf;
plot(w_uniform);
plot(w_sinc,'r.');
title('Uniform weighted (blue solid) and sinc weighted (red dotted) sensor weightings');


Hilbert transform with strange behavior I am taking the Hilbert transform of a signal in Matlab and I don't understand part of the result.

Below are the complete signal $z(t)$ and phase $\phi(t)$, obtained with

phi=unwrap(angle(hilbert(z)))

I noticed a jump around point 400, so I zoomed this region. This is shown below.

Although the signal in this region is very well behaved, the phase has an almost discontinuity. Why?

It is not due to the unwrap command, as we can see below

The data can be found here

A Hilbert Transformer is an acausal linear time-invariant (LTI) filter.  Being LTI it has all of the things an LTI filter has: the output is the convolution of the input against an impulse response.  And it has a frequency response with magnitude and phase.

But being acausal, that means that the impulse response begins to respond before the impulse hits the input (i.e. it can see future input).  That, of course, is impossible in a real-time system so the impulse response of the Hilbert Transformer is delayed to make it causal.  That delay causes a linear trend to the phase component.

The magnitude frequency response of the Hilbert Transformer is 1 (as linear gain) or 0 dB for virtually all frequencies (DC and Nyquist are the exceptions). But without the delay, a Hilbert Transformer shifts the phase for all positive frequencies by -90° (or $-\frac{\pi}2$ radians) and all negative frequencies by +90°.  That puts a phase discontinuity at 0 Hz (DC).  Likewise, you will see the opposite phase discontinuity at Nyquist.


general purpose microprocessor vs DSP(digital signal processor)? I know that DSPs are specialized processors but i want to know what will happen  if DSPs are replaced by general purpose processor in their applications? I know one application,mobile phone,will that work normally if DSP of mobile is replaced by a normal(general purpose microprocessor)

On a basic level, in my opinion, a DSP chip must, as a minimum requirement, be able to optimally perform the convolution summation used to compute the output of an FIR filter.
$$ y[n] = \sum\limits_{k=0}^{K-1} h[k] \, x[n-k] $$
$y[n]$ = FIR filter output
$x[n]$ = FIR filter input
$h[k]$ = impulse response of FIR filter (also FIR coefficients)
To do this in $K$ instructions, one must be able to in one instruction:

Do a multiplication and addition, implemented as a single multiply-accumulate (MAC) instruction.

Fetch both the delayed sample, $x[n-k]$ and the FIR coefficient $h[k]$ for the next MAC.  To do both, this requires at least two data memory spaces: one for coefficient array $h[k]$, one for the signal array $x[n]$.  Both spaces can be accessed in parallel (simultaneously).  This is called a "Harvard Architecture".

Be able to saturate the results to defined maximum positive or negative values if the result exceeds that range.

Besides being able to saturate to the "rails", the Arithmetic Logic Unit (ALU) in the DSP should be able to round to the nearest LSB of the high-order bits (that will be output or saved to memory) based on the low-order bits of the accumulator.

To address the delayed samples in memory, the ability to do modulo arithmetic on the pointer to $x[n-k]$ is necessary.  This implements a "circular queue" or "circular buffer" or "first-in-first-out" (FIFO) buffer which, in the audio DSP world, is just called a "delay line".  The DSP programmer should not have to worry about when a pointer increments beyond the boundary of a circular queue requiring the pointer to be adjusted or wrapped around.  The DSP chip should make that modulo adjustment to the pointer automatically.

In addition, for the FFT, the ability to do bit reversal or bit-reversed addressing is necessary along with some other handy instructions that make for efficient coding of FFT butterflies.

Now, in order to be a "DSP chip", it really should be able to do all of these operations without explicit instructions in the inner loops.  They build it into the hardware of the ALU the means to perform those operations without explicit effort from the programmer (nor extra machine instruction cycles).  Nothing prevents you from using a normal CPU to do all of these operations, even in a "high-level" language like C/C++ .  And the general purpose CPU chip might be running at 5 GHz while the DSP chip is running at 500 MHz, so the 10x instruction speed might more than make up for the cost of all those extra instructions that the algorithm programmer must code.


How do optical anti-aliasing filters work from a frequency domain perspective? To prevent aliasing caused by the finite number of pixels on a sensor, a blurring filter is commonly used. How does that work from a frequency domain perspective? What is the transfer function of such a filter?

In Karel Fliegel, Modeling and Measurement of Image Sensor Characteristics, Radioengineering, vol. 13, no. 4, December 2004, he gives the optical transfer functions (OTFs) for different detector photosensitive area shapes. OTF is the Fourier transform of the spatial domain impulse response. The spatial domain impulse response, also called the point spread function (PSF), is a 2-d function that is constant-valued (assuming uniform sensitivity) inside the detector photosensitive area and zero outside it.

If there is an optical low-pass (blurring) filter in front of the sensor or implemented by shaking the sensor, another PSF describes it, and the total PSF is the convolution of the PSF of the sensor element and the PSF of the optical low-pass filter. The total OTF is the product of the two OTFs. Birefringenge-based optical filters (for example patent US 6937283 B1) create shifted copies of the image so their PSFs are sums of shifted and possibly scaled Dirac pulses. For example, a square arrangement of four Dirac pulses is equivalent to a cascade of a horizontal and a vertical comb filter.

An optical anti-aliasing filter also attenuates in-band frequencies that would not alias. Using a software sharpening filter the amplitudes of those frequencies can be restored. This also amplifies aliased frequencies, but because their attenuation in the combined OTF of the anti-aliasing filter and the detector is stronger, they retain some of that attenuation.

Here is an example analysis of a horizontal-vertical-separable total PSF that includes 1) a square photosensitive area, 2) a birefringent optical low-pass filter that has a PSF that consists of a square arrangement of Dirac pulses separated horizontally and vertically by 0.4284281154 times pixel width, and 3) a digital post-processing filter (operating on the pixel data), which can be separated to a cascade of a horizontal and a vertical filter, both with a PSF [-0.01536945896, 1.030738917, -0.01536945896]. A square grid of pixels with no gaps between pixels (100 % coverage of the imaging plane by photosensitive areas) is assumed. The analysis is done only along one of the separable dimensions for simplicity, but it would be better to do a 2-d analysis.

Figure 1. Red: combined one-dimensional PSF of the optical anti-aliasing filter, the square photosensitive area, and the digital sharpening filter. Turquoise: the ideal sensor PSF for a square grid (sinc function). Black: PSF of the square photosensitive area. Horizontal axis: displacement in units of one pixel width.

Figure 2. One-dimensional OTFs (Fourier transforms) of PSFs of Fig. 1, using the same color code. Horizontal axis: spatial frequency in units of radians per pixel width. There are some plotting artifacts that make parts of the red curve invisible.

Considering sinc (Fig. 1) as the ideal one-dimensional PSF, the mean square error for imaging spatial white noise is reduced by 40 % by using the optical anti-aliasing filter and the digital sharpening filter. The sum of squares error can be calculated as the integral over the square of the difference between sinc and the actual PSF, or by Parseval's theorem equivalently in the frequency domain using the differences of the OTFs (Fig. 2). I optimized the coefficients and the amount of separation in the optical anti-aliasing filter so that the mean square error is minimized.

Without the optical anti-aliasing filter, a digital post-processing filter (not a sharpening filter this time) with optimized coefficients [0.06765983302, 0.8646803339, 0.06765983302] only reduces the mean-square error by 11 % compared to a reduction of 40 % from the combination of an optical anti-aliasing filter and a post-processing sharpening filter, so it is apparent that optical anti-aliasing is useful.


Bandwidth of an energy signal I'm studying for my midterm , and I'm struggling with finding bandwidths.

I have this energy spectral density : $\Psi(f)= { 1  \over \sqrt {  1 +  ({f \over B }) ^2} }$ ,  and I need to find  the $3$ $dB$ band and the band with $90%$ % power.

About the $90%$ % case, do I need to find the total energy and then  $\int_0^W \Psi(f) = 0.9 \times  Total Energy$?

One think that doesn't make sense for me is that I'm dealing with an energy signal but I need to find something related to power.

If I had a power signal and wanted to find the 3 dB bandwidth, would I need to find the average power and then solve $\int_0^W S_x(f) = 0.5 \times TotalPower$ ?

Thanks in advance, and sorry for the bad english.

I'm using the following assumption: the magnitude of the spectrum is given by the formula

$$|X(f)|=\frac{1}{\sqrt{1+(f/B)^2}}$$

Then the energy spectral density is

$$|X(f)|^2=\frac{1}{1+(f/B)^2}\tag{1}$$

The 3 dB bandwidth is simply the frequency for which (1) becomes $0.5$, which is $W_{3dB}=B$. The bandwidth according to the 90% energy/power definition is

$$\int_{-W}^W|X(f)|^2df=0.9\cdot\int_{-\infty}^{\infty}|X(f)|^2df\tag{2}$$

Equation (2) is the reason why I assume the definition as given in (1), because otherwise the integral on the right-hand side of (2) does not converge, so the 90% energy/power definition would not make sense. With the definition (1) the integrals in (2) are easily solved using the $\arctan$ function, from which $W=const\cdot B$ follows, where $const=6.31$ (if I'm not mistaken), but I'm sure you can do the calculation yourself.


Sketch the following Signal: $x(t)=3\sin(\frac{\pi t}{4})$ I feel kind of overwhelmed by it, and would like someone to walk me through the steps of sketching the signal.

From your physics knowledge you get that sine waves are usually equated, in the general case, as

$$y(t) = A * \sin(\omega t)$$

where $\omega = 2 \pi f$ and the periodicity is $T=1/f$

Comparing the general equation with your particular case, you have: $$2*\pi*f = \pi / 4$$
meaning that $$f=1/8 Hz$$
or $$T=1/f=8s$$
The amplitude is $A=3$

If you want to plot it in matlab, for example, just write

t = 0:0.01:16;
x = 3*sin(pi*t/4);
plot(x,t)

As a final comment, this is really a basic question and you should carefully read your textbook if you want to progress to more advanced topics.


Can every type of linear filter be modelled by a convolution? I have an input time series going through a filter that creates another time series as output.

If I assume in first approximation that my filter is linear, does it necessarily mean that I can model this filter by a convolution?

In other words is it only certain types of linear filters that can be modelled by a convolution? (LT1?...).

I have never seen any clear statement on this in any publication. It is as obscure and vague as definitions of time invariance, ergodicity, causality etc in filtering textbooks.

My question comes from wondering if there is any point trying to attempt an already risky deconvolution if I am not even sure that my linear filter is necessarily akin to a convolution.

Why classes of linear filters are not clearly explained in publications (that I know of) is a mystery to me, as it should perhaps be the first question that arises when trying to model real world phenomenons.

No. It's only LTI (Linear and Time-Invariant) systems that can be modeled with convolution through a unique single impulse response.

For example the systems

$$ y(t) = g(t) x(t) $$ or
$$ y[n] = \sum_{k=0}^{k < n} x[n-k] $$

are both linear but not time-invariant and their output $y[n]$ cannot be computed with the convolution operation ( $\star$ denoting convolution)

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

as such an impulse response $h[n]$ does not exist.


Constructing frequency filter from a 2d spatial filter Given a 2d spatial kernel $h$, I want to construct the corresponding frequency filter $H$. In theory, $H$ is simply obtained by calculating DFT of $h$ using a FFT algorithm. But in practice, there are some considerations to be considered. First, if $H$ is used for filtering, zero-padding is required. Second, after zero padding, the original center of $h$ (i.e. center of $h$ before zero padding) must circularly-shift to the upper-left corner.
I know that zero padding is used to overcome wraparound errors. And I have observed, by experiment, shifting the original center to the upper-left corner prevents shifting the filtered image. But my question is why this is the case? That is, why not shifting the original center, leads to shifting the final image; and how shifting the original center resolves the problem?
Thanks

I found it! Placing the original center at the upper-left corner, by circularly-shifting $h$, makes $h$ an even function. From the properties of Fourier Transform, we know the DFT of a real and even function is real. So, $H$ becomes real and its phase angle is zero. Hence, the phase of the filtered image is as the phase of the original image.
By the way, shifting $h$ does not affect the frequency spectrum of $H$.


Is it possible to design a band-pass filter with fewer taps than a series HPF + LPF? I am digitally filtering a signal that is streaming live. The goal is to reduce the filter delay as much as possible while maintaining the quality of some measurements based on the filtered signal. For legacy reasons, the system was built with a low-pass and high-pass FIR filter in series - each filter has the same number of taps (kernel size), $n$.

I realized though, we are no longer constrained to these and could implement a single FIR band-pass filter. If I can reduce number of filter taps involved from $2n$ that will reduce the delay. My question is: is it possible to design this band-pass filter with fewer taps (a smaller kernel) than the combined taps (kernels) of the LPF/HPFs while achieving the same performance?

My initial thought is FIR filtering is a convolution operation, and thus linear + associative. That means you can construct the band pass kernel by convolving the LPF + HPF kernels. I'm unsure how to interpret this in terms of taps/kernel size because the convolution in some sense assumes an infinite kernel.

Thanks!

Assuming that your legacy FIR linear phase filters were designed using windowing method, then yes, you can combine those two FIR linear phase filters into a single FIR linear phase filter and it will be shorter than the the serial cascade.

However, you won't be convolving the individual impulse responses to get it. Rather, you should use a multi-band windowed FIR filter design technique. An example of how to do it can be found on Alan Oppenheim's DTSP book on chapter 7 filter design section.

Specifically he makes use of Kaiser window design approach acording to which the filter (window) length is empricially dependent on the transition band width and peak approximation error.

So, given the same transition band width and error ripples of a HPF -> LPF cacade, whose lengths are M samples each, you would be using a single multiband filter of length M with the same transition band and error ripple.

Note that as MattL also indicates, you may also chose different filter types to achieve similar respones, but that depends on what you want to achieve finally.


different moving average filters differences Suppose we want to use a moving average filter (3-point moving average) for $X[n]$ (where $Y[n]$ is output). What are the differences among

$Y[n]=\frac13X[n-1]+\frac13X[n]+\frac13X[n+1]$,
$Y[n]=\frac13X[n]+\frac13X[n+1]+\frac13X[n+2]$, and
$Y[n]=\frac13X[n]+\frac13X[n-1]+\frac13X[n-2]$ ?

Hi: All of those moving averages require future values of $X$ which in DSP terminology means that they are non causal filters and physically not realizable. Also, in general you should not have $X[n]$ on both sides of the equation because simplification is then possible and this means that it's no longer an equally weighted filter.  So, if you fix that and move the index back one so that you have $Y[n] = 1/3 *(X[n-1] + X[n-2] + X[n-3])$ then that is a causal filter that is realizable.


what are the vector formats similar to hash codes to store user's voice for future voice verification? I have wav file of ~2-10sec of user's voice, taken from a conversation.
Maybe it's SIMILAR TO a hash code.
the general voice verification should look like:

wav file -> vector
new wav file ->vector
compare 2 vectors with cosine similarity

which is the convenient way to store wav file to be able to verify it later on with another wav voice file?

Since you are not trying to verify if two "messages" are identical in a word-to-word sense, it would not be practical to go down the hash-function way.

This sounds more like a "pattern recognition" problem.

The typical approach in this kind of problems is to obtain a set of features from the signal and then, given the features of some query signal, see if they match. The matching there is more in the sense of Euclidean Distance between vectors of features. Of course, other metrics of distance can be used too. This is captured in a straightforward way in the case of the Nearest Neighbour algorithm.

Usually, in the case of speech signals and voice recognition, the Mel Frequency Cepstral Coefficients are used as this vector of features you are after.

You might also find this response useful, as well as this repository.

Hope this helps.


Why is there a difference between MATLAB's FFT of a zero-padded signal compared to the non-padded case I created a data signal and calculated the spectrum via the fft-command in MATLAB. The data-signal itself is not periodic but it is created in a way that it can be treated as one period of a periodic signal. The last sample smoothly transitions to the first one. I observe that the absolute of the spectrum strongly depends on the zero-padding before calculating the fft. If I do zero-pad the signal such that its length is a power of two, the spectrum differs severely from the one where I kept the signal unchanged. In the following pictures you can see the results.

It can be clearly seent that there is a kind of gap at the center of the channels without zero-padding. I can't explain this. I would expect the shape of both transforms to be the same. Why is that not the case here, does anyone have an idea?

If you write

FFT(data)

you basically sample a continous "spectral curve" at length(data) many points. If you do zero padding

FFT(data,len)

where length(data)<len you effectively sample the curve at more points (len) that are more closely spaced. But it is still is the same curve. It was just sampled at different frequencies.

If this is really the only difference you made, then you're probably witnessing a "sampling issue". By that I mean the curves appear to be different because you've sampled it at different points and there is some oscillations going on. Chances are that you havn't applied any windowing which would smooth out the curve a bit and make this effect (otherwise known as spectral leakage) less likely.

You could try one of the following things:

do even more zero padding for an increased sample rate (samples/Hz) and check out the shape
apply a window function to smooth out the curve before the FFT
use PWELCH

PWELCH does more than simply an FFT. It includes windowing and averaging of overlapping but smaller windows. This is probably what you should be doing given your "non-periodic" hint.


Linear Predictive Coding for general signals I have a signal that's monotonic and roughly linear and have been looking at using Linear Predictive Coding to encode information and compress my signal. I guess my first general question is if this is even possible for general signals? I also see LPC mentioned in the context of speech so I'm wondering if LPC is particularly effective for one type of signal and would be really bad for another type of signal. Here's how I understand the process so far:

We have an input signal $x(n)$ that we look to model as $$x(n) = \sum_{i=1}^{p} a_i x(n-i) + e(n)$$

where the $a_i$ are our prediction coefficients and $e(n)$ is the error. We can find the error $e(n)$ by filtering the inverse transfer function $\frac{1}{H(z)}$ with our signal $x(n)$. Thus, if we can closely model the error with a small number of a parameters, our signal has been reduced to our prediction coefficients and these small parameters.

I see why LPC is often used for speech compression - because the error $e(n)$ is roughly white and so we can roughly model it by its variance. However, my intuition tells me that it should be much easier to compress something like perfectly linear signal than a highly nonstationary signal like a vocal sample. So I guess my questions are (provided I'm understanding LPC correctly):

Do people use LPC for things other than speech, and if so, how do they go about modeling $e(n)$ to achieve compression?
If my intuition is correct about it being easier to compress a linear signal using LPC, how would one do it in this case? I'm confused because since a perfectly linear signal is exactly predictable, it should produce $0$ for all $e(n)$. However, filtering a linear signal with prediction coefficients $[2,-1]$ produces a non-zero $e(n)$.

There may be a misunderstanding about LPC here. For one, it is mainly used for speech because it was designed to fit the model of the human vocal system. The figure below is taken from Purdue's VISE lab in one of their speech processing exercises. The figure shows speech modeled as a filter that is excited by either an impulse train or white noise. LPC works well for speech because it gives you a concise way of describing the filter part. You also need to estimate the driving frequency of the pulse train, detect if a time segment is voiced or unvoiced, and control for the power of each time segment.

In the frequency domain, the a speech signal will look like

and all the LPC is expected to do is to fit the spectral envelope.

  I see why LPC is often used for speech compression - because the error e(n) is roughly white and so we can roughly model it by its variance.

The significance of a white, particularly additive white Gaussian, error signal is that it means that you have used all of the information. If your error is highly structured, it means that you could do a better job at modeling the system.

  Do people use LPC for things other than speech, and if so, how do they go about modeling e(n) to achieve compression?

I'm not aware of it being used for anything apart from speech. As stated above, the vocal system is being modeled by LPC, not the error. The vocal system is modeled using the source-filter model, and LPC is a good way of describing the filter part of the model.

I'm not sure about your second question, but I hope the rest of this help. Speech encoding with LPC is a really interesting, particularly because it works so well. If you want to do some matlab experiments, you could try Lab 9b on Purdue's VISE lab.


Algorithms for finding fundamental frequency based on ACF result I am developing a software for fundamental frequency tracking. For this purpose, I have designed a function which calculates autocorrelation over the signal and a second function which, based on autocorrelation result, finds the peak corresponding to the fundamental frequency.

However the code I have written is not very reliable. What I want to know is: which algorithms are reliable for fundamental frequency peak search?

EDIT

My autocorrelation plot looks as follows so far, when done in MATLAB:

however, I implemented it in C, because my software is in C, as follows:

 91     for(;i < *length;i++) {
 92         for(j = i;j < *length + i;j++) {
 93             if(j > i) {
 94                 sum += 0;
 95             } else {
 96                 sum += samples[i] * samples[j];
 97             }
 98         }
 99
100         result[i] = sum;
101         sum = 0.0f;
102     }

with i starting in 0. I also believe this code to be OK.

Approximating the above plot in MATLAB, I see this, where the marked peak would be my fundamental frequency one (the note I'm testing with is an E2)

as we can see from this plot

the signal I'm dealing with is quasiperiodic and the difference between those peaks actually result in the fundamental frequency for E2 note: the distance between the two is 389 - 193 == 196 (values obtained from drag & drop of MATLAB plot), which, divided by my sampling rate of 16 kHz, results in 0.01225, that would be my wave period. Inverting this value I obtain 81.6326530612 Hz, which is very near to my expected 82.41 Hz for E2.

However I'm having a real bad time trying to locate these two peaks in code, since I'm not sure what to compare to be sure the peaks I found are actually these two peaks.

im assuming your signals are audio and harmonic and you are not doing multiple F0 estimation where you have multiple tones played simultaneously. in that case there are a few easier methods. the normalisation of short time autocorrelation function. r(lag) = (1/N)*sum(x(n)*x(n+lag)
N is length of x(n)
F0 here is obtained from the inverse of the lag that corresponds to the maximum of r(lag) within a predefined range (for example your range could be from 20 to 250, hence at lag 196 you have a maximum, again play around with this range). to avoid detecting an integer multiple of the period, short lags are better than long ones.
here is a C function for this, probably worth playing around with.

void autocorr(
    float *r,       /* (output) autocorrelation vector */
    const float *x, /* (input) data vector */
    int N,          /* length of data vector */
    int order       /* largest lag for calculated
                      autocorrelations */
)
{
    int     lag, n;
    float   sum;

    for (lag = 0; lag <= order; lag++) {
        sum = 0;
        for (n = 0; n < N - lag; n++) {
            sum += x[n] * x[n+lag];
        }
        r[lag] = sum/N;
    }
}

An alternative method for detecting peaks, is center clipping. after your autocorr function you could code something like:

y(n) = a(n) - C,  if a(n) >= C

y(n) = a(n) + C, if a(n) <= -C

y(n) = 0 otherwise.

C is some percentage of the max amplitude of the autocorr vector.

all i have showed you above is the very basic stuff, since F0 analysis can be very complex (probably why people have not answered). you can find more in speech codec books or music transcription books.


Create 16 bit sine and write to file This seems a simple task, yet I am making some stupid mistakes I am not able to find. I create a sinus wave (integer) and then convert it to 16 bit char, little endian. Then I save to file and open it in Audacity (see image below). The first channel in the image is created by Audacity, the second by me. Something is really wrong there. I think it must be in the bit shifting, because if I use 8 Bit, it works. So the extra byte causes trouble but I don't see, how.

#include <bitset>
#include <iostream>
#include <chrono>
#include <math.h>
#include <fstream>
#include <iomanip>

#define PI 3.14159265359

void
test_sine16Bit(){

    int samplesize          = 16;
    unsigned int samplerate = 44100;
    int duration            = 3; // seconds
    int channels            = 1;
    //long samples            = 4096;
    long samples            = samplerate * duration; // per channel
    long total_samples      = samples * channels; // all channels

    // Create test input buffer
    char input[total_samples * samplesize/8];

    // Initialize test buffer with zeros
    for (int i=0;i<total_samples * samplesize/8;++i){
        input[i]=0;
    }

    // Create test input
    int total_index=0;
    char graph[20];
    for (long i=0;i<samples;++i){
        // Channel 1: 440 Hz sinus
        int sinus    = .8*0x8000*sin(2*PI*440*i/samplerate); //  amplitude = 0.8 * max range; max range = 0x8000 = 32768 ( max value for 16 Bit signed int )

        // Little Endian: BUG HERE?
        input[total_index] =  sinus & 0xFF;             // LoByte
        //input[total_index+1] = (sinus & 0xFF00) >> 8; // HiByte
        input[total_index+1] = (sinus >> 8) & 0xFF;     // HiByte
        total_index += 2;

    }

    // Export raw audio data to file
    std::ofstream fout("sine_char16bit.raw");
    if (fout.is_open()){
        long idx = 0;
        for(long i = 0; i < samples; ++i){
            fout << input[idx] << input[idx+1]; idx+=2;
        }
        fout.close();
    }
}

EDIT:

Found the solution: You need to open the file in binary mode and not write via <<, but with put or write.

std::ofstream fout("sine_char16bit.raw", std::ios::out | std::ios::binary);
if (fout.is_open()){
    fout.write(input, total_samples * samplesize/8);
}
fout.close();

Create low-frequency triangular or square waveform. Start from constant value, written into the wave. Generate one in audacity and save. Compare. Do not generate sine or anithing more compled right away. Start building complexity from simple components. Be smart, use simplicity and evolutionary approach.


Database into training and testing images I am trying to compute the accuracy of the HMAX model. I am using the Face category (containing 435 images) from the Caltech101 database. I split it into x training and y testing. At each time, when x increases, the accuracy also increases. Furthermore, I heard that the number of training should be equal to 80% by comparing it to the tests. So when I split my data into 348 positive training and the rest for positive testing, I got an accuracy that it is smaller than the other smaller splits (when x<348)!!.

Even when I took x=300, I also got a smaller accuracy!!

By the way, I also used the background category and I split it into 50 negative training and 50 negative testing.

Please why I got a smaller accuracy? Please I need your help.

Some of the possible reasons would be:

Object recognition accuracy is not really well-defined. Instead accuracies are for particular datasets.
What you are working with is called cross-validation, where you split your own dataset into testing and training datasets and then use the trained model to make predictions on the test dataset.
Depending on how you do the cross-validation and the final model, the accuracy tests may vary. For example, if (even by a random process) your training set contains well-conforming images and the test dataset contains images which are not even closely conforming to your model, then the accuracy would be low.

As for a solution to your dilemma, you should probably consider looking at the dataset and sampling in such a way that this situation can be avoided. Alternatively you could use a more generalizable algorithm and check for overfitting of the data via cross-validation again.


Inverse Fast Fourier Transform Hard Way I'm reading "Digital Signal Processing", written by Michael Weeks. In section Representing a Digital Signal as a Sum of Sinusoids, there is this Matlab code that I cannot understand. (code given below)

What I don't get is what does the term inside cos argument 1/Xsize mean (in part %Do the sinusoids as discrete points only, I guess it is fundamental frequency, but how can we conclude it). And why is the whole function divided by Xsize.

Same thing with smooth0, smooth1, smooth2... I also don't get this smoothing process.

Thank you for any explanations, or links to helpful stuff.

%sum_of_sins.m
%Show a random signal as a sum of sins
%

%Make our x signal
x=round(rand(1,8)*10);
Xsize=length(x);

%Get FFT of x
X=fft(x);
Xmag=abs(X);
Xphase=angle(X);

%Show the freq-domain info as sum of sinusoids
%Find the IFFT (the hard way) part 1
n=0:Xsize-1;
%Do the sinusoids as discrete points only
s0=Xmag(1)*cos(2*pi*0*n/Xsize+Xphase(1))/Xsize;
s1=Xmag(2)*cos(2*pi*1*n/Xsize+Xphase(2))/Xsize;
s2=Xmag(3)*cos(2*pi*2*n/Xsize+Xphase(3))/Xsize;
s3=Xmag(4)*cos(2*pi*3*n/Xsize+Xphase(4))/Xsize;
s4=Xmag(5)*cos(2*pi*4*n/Xsize+Xphase(5))/Xsize;
s5=Xmag(6)*cos(2*pi*5*n/Xsize+Xphase(6))/Xsize;
s6=Xmag(7)*cos(2*pi*6*n/Xsize+Xphase(7))/Xsize;
s7=Xmag(8)*cos(2*pi*7*n/Xsize+Xphase(8))/Xsize;
%Redo the sinusoids as smooth curves
t=0:0.05:Xsize-1;
smooth0=Xmag(1)*cos(2*pi*0*t/Xsize+Xphase(1))/Xsize;
smooth1=Xmag(2)*cos(2*pi*1*t/Xsize+Xphase(2))/Xsize;
smooth2=Xmag(3)*cos(2*pi*2*t/Xsize+Xphase(3))/Xsize;
smooth3=Xmag(4)*cos(2*pi*3*t/Xsize+Xphase(4))/Xsize;
smooth4=Xmag(5)*cos(2*pi*4*t/Xsize+Xphase(5))/Xsize;
smooth5=Xmag(6)*cos(2*pi*5*t/Xsize+Xphase(6))/Xsize;
smooth6=Xmag(7)*cos(2*pi*6*t/Xsize+Xphase(7))/Xsize;
smooth7=Xmag(8)*cos(2*pi*7*t/Xsize+Xphase(8))/Xsize;
%Find the IFFT (the hard way) part 2
my_sum_of_sins=s0+s1+s2+s3+s4+s5+s6+s7;
smooth_sum=smooth0+smooth1+smooth2+smooth3+smooth4+smooth5+smooth6+smooth7;

%Show both discrete points and smooth curves together
xaxis1=(0:length(smooth0)-1)/20; %for 8 points
xaxis2=(0:length(s0)-1);

figure(1);
subplot(4,1,1); plot(xaxis1, smooth0, 'g', xaxis2, s0,'r*');
subplot(4,1,2); plot(xaxis1, smooth1, 'g', xaxis2, s1,'r*');
subplot(4,1,3); plot(xaxis1, smooth2, 'g', xaxis2, s2,'r*');
subplot(4,1,4); plot(xaxis1, smooth3, 'g', xaxis2, s3,'r*');

figure(2);
subplot(4,1,1); plot(xaxis1, smooth4, 'g', xaxis2, s4, 'r*');
subplot(4,1,2); plot(xaxis1, smooth5, 'g', xaxis2, s5, 'r*');
subplot(4,1,3); plot(xaxis1, smooth6, 'g', xaxis2, s6, 'r*');
subplot(4,1,4); plot(xaxis1, smooth7, 'g', xaxis2, s7, 'r*');

Effectively Xsize is allows the sampling frequency to be set.

If you just used 2*pi*0*n+Xphase(1) as the argument of $\cos(\cdot)$ or $\sin(\cdot)$ then you'd just get $\cos({\tt Xphase(1)})$ or $\sin({\tt Xphase(1)})$ for all terms... because then the 2*pi*0*n terms would be an integer multiple of 2*pi.

EDIT
With respect to the Xsize terms at the end:

  s0=Xmag(1)*cos(2*pi*0*n/Xsize+Xphase(1))/Xsize;

Check out the definition of the inverse DFT, particularly equation 2:
$$
x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k e^{i2\pi k n / N}
$$
for your code, N = XSize.


Variations in MATLAB plot each time the program is executed Every time I run the program my PLOT in MATLAB gets changed to some extent (some times to larger extent too).
I understood that the channel is random and this might happen due to random distributions used to model the channel. Is there any way to solve this issue, so that every time I get the same values and plot.

The plots are changing i.e. you are changing inputs to the model (be it random input sequence , random channel matrix , random gaussian noise ) .. there are lot of random combinations here.. which will definitely change the plot.. you cannot expect to have same plot everytime when you are dealing with random numbers..
You can control with using same seed in random number generator.

Use the same seed while dealing with random numbers

Quick google search on this gave me this which you can use

https://www.mathworks.com/help/matlab/math/generate-random-numbers-that-are-repeatable.html


Draw a frequency response graph from a simple equation How is it possible to draw an audio frequency response graph (x-axis : freq from 20 to 20 000hz, y-axis : dB)  from a simple equation like

y[n] = x[n] - x[n-486]

(Example: with this equation I get a nice comb filter that removes frequency 197,5Hz and all its harmonics, for an audio file with sampling rate 96 Khz)

How to turn such an equation into a graphical interpretation like this
http://downloadfreesamples.com/wp-content/uploads/2013/08/EQ-Four-CIS-DSP_3.jpg ?
Then I could see how deep / steep the filtering is.

Thank you

(PS : I continue to read references and online books about filter design that were given to me in previous questions, thanks for these references by the way.)

The tool of choice here is the z-transform, which is applied to your equation. It transforms sequences of complex numbers to functions in the complex plane. Its most useful properties are that it takes delays to powers of the argument z, and for z = exp(i omega) we get the fourier transform of the transformed sequence. Look up the details, I'll just show you how it works here.

y[n] = x[n] - x[n-486]  ---> Y(z) = X(z) - X(z)*z^(-486)

The equation in z-domain factors on the right hand side to give Y(z) = X(z) ( 1 - z^(-486) ).

Since Y(z) is the output and X(z) is the input, for linear systems we define the transfer function to be Y(z)/X(z) = H(z) and the above equation becomes

H(z) = Y(z)/X(z) = 1 - z^(-486)

Now if we substitute z = exp(i omega) we get the frequency response of your system, which is

H_f(omega) = 1 - exp(i omega)^(-486)

and with the properties of exp() it simplifies to

H_f(omega) = 1 - exp(- 486 i omega)

If you need actual frequency f instead angular frequency omega, just substitute omega = 2 pi f/fs, with sampling rate fs.

Finally, you'll probably want to know just the magnitude of the response, so look at abs(H_f).
In this specific case it's very simple to calculate because abs(H_f) = sqrt( H_f * conj(H_f) ) and conj(H_f) = 1 - exp(+ 486 i omega).

As requested, some code. Here's what you can do in matlab/octave to get the magnitude frequency response with a logarithmic frequency axis and magnitude axis in dB.

% Create a logarithmically spaced frequency axis with range 20Hz to 20kHz

fAxis = logspace( log10( 20 ) , log10( 20000) , 10000);

% calculate the complex frequency response at these frequencies

% Set SampleRate before calling this, like SampleRate = 44100;

Hf = 1 - exp(1i*486*2*pi*fAxis/SampleRate);

% plot the magnitude spectrum

semilogx(fAxis,10*log10(Hf.*conj(Hf)));


Proving that the IDTFT is the inverse of the DTFT? The DTFT is given by:

$$X(e^{j\omega}) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}$$

The IDTFT is given by:

$$x[n]=\frac{1}{2\pi}\int_{0}^{2\pi}X(e^{j\omega})e^{j\omega n}d\omega$$

I have been able to show by substitution of the DTFT into the IDTFT that the transform and a subsequent inverse transform return $x[n]$:

$$\begin{align}
x[n]&=\frac{1}{2\pi}\int_{0}^{2\pi}X(e^{j\omega})e^{j\omega n}d\omega\\
&=\frac{1}{2\pi}\int_{0}^{2\pi} \left( \sum_{k=-\infty}^{\infty}x[k]e^{-j\omega k} \right)e^{j\omega n}d\omega\\
\end{align}$$

Swap the order of integration and summation:

$$x[n]=\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}\int_{0}^{2\pi}x[k]e^{j\omega (n-k)}d\omega$$

Argue that given $e^{j\omega (n-k)}$ is an odd function, it will only evaluate to anything other than 0 when $k=n$:

$$\begin{align}
x[n]&=\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}\int_{0}^{2\pi}x[k]e^{j\omega (n-k)} d\omega \ \delta[n-k]\\
&=\frac{1}{2\pi}\int_{0}^{2\pi}x[n]d\omega \\
&=\frac{2\pi}{2\pi}x[n]\\
\end{align}$$

However, I have been unable to show the dual case: that the inverse transform (IDTFT) substituted into the forward transform (DTFT) gives $X(e^{j\omega})$. How can we show this?


A simple (hanning-windowed) FFT reveals four peaks that stand out far beyond the noise. This appears to be what you are looking for.


Adding attenuations in dB I am writing a synthesizer. One part of the synthesizer specification says that different sources can specify an attenuation to be applied to a signal, and that these attenuations should be summed.
The attenuations are specified in cB and are given as positive values. An attenuation of 0cB means no attenuation.
The naive approach, which I've tried, would be:
$$
x_{total} = 20 \log_{10}\left( \sum 10^{\frac{x}{20}}\right)
$$
for each $x$ as the attenuation (converted to dB) being summed.
This has problems, since if I add two 0cB attenuations (i.e. no attenuation, twice), I get a total attenuation of 60.2cB. This is obviously not what I want.
How can I add my attenuations to get a sensible result out the other side? Thank you.

Part of the advantages of using dB to express gain or attenuation is that they're inherently logarithmic.
So if you're speaking in dB, when you add a 10dB attenuation to another 10dB of attenuation, what you're saying is that you're multiplying the signal by $\sqrt{0.1}$, and then by $\sqrt{0.1}$ again.  The result is 20dB of attenuation -- or multiplying the signal by $\sqrt{0.01} = 0.1$.
So if your ultimate gain is $A$, as a normal number, and you have attenuations $a_1, a_2, \cdots$ in dB, then the actual $A$ by which you multiply your signal is $$A = 10^{-\frac{1}{20}\sum a_n}$$


How to determine where to sample for demodulation of BPSK signals? I have a simple BPSK demodulator. Very simply, the signal comes in and is split into two branches, one for I and one for Q.

The I branch is mixed with a sin wave of the carrier, and the Q branch mixed with the cosine wave of the carrier.
Then, each of the outputs is convolved with a matched filter, (in this case a root-raised-cosine, which initially pulse shaped my transmitted bits).
Now, if I combine both outputs and take the magnitude of both of them (I^2 + Q^2), I get the envelope.
From here, I look for peaks, and choose their indicies.
I then use those same indicies that I got from (4), and sample my I and Q at those same indicies. Now I have complex softbit values, and am able to correct for any phase or frequency offsets, and can make hard decisions from there. (I have the luxury of post processing btw).

This technique worked well with high to medium SNR. The problem is that with lower SNR (and maybe multipath), I do not know how to perform bullet (4). The problem is that instead of one peak in the envelopes corresponding to one bit, I now have multiple peaks, and 'peak picking' wont work because there are so many to choose from... how is this problem mitigated? This seems to be a bottleneck because if I cannot tell where my 'bit' is, I cannot move forward. Any thoughts?

Any help is appreciated, math is good, although I would ask the answers be more pesudo-code /concept oriented in advance. Thanks!

Edit: Thanks to feedback from Dilip, I should also add that I can clearly see 'energy' in the envelope when there is a signal and there isnt a signal. (ie, before my signal arrived), so this leads me to believe that I should be able to sync/demod and that it isnt an SNR issue solely...

I have done something similar to this in MATLAB. In my cause, I used an Early/Late Gate clock recovery method to get an estimate of the offset between transmit and receive symbol timing. This method uses 3 samples per symbol - one at the optimal sample time, one that is 1 sample delayed and 1 that is one sample advanced. This works well for on/off keying signals, since there is no negative portion of the signal, but the methods below I think will work better for bipolar signals.

Another method is the Muller and Mueller Algorithm, which uses 1 sample per symbol. This uses the equation $e_n = (\widehat{y_n}y_{n-1}) - (y_n\widehat{y}_{n-1})$, where $e_n$ is the error , $y_{n-1}$ and $y_{n}$ are the previous and current symbol samples, and $\widehat{y_n}$ and $\widehat{y}_{n-1}$ are the decision outputs (-1 or +1 in your case). M&M requires that the carrier is recovered first.

A third algorithm, which is very similar to the one above, is the Gardner algorithm. In this algorithm, the error signal is $e_n = (y_n - y_{n-T})y_{n-T/2}$, where $y_{n-T}$ is the sample one symbol in the past, and $y_{n-T/2}$ is half a symbol in the past. This requires 2 samples per symbol, but is not sensitive to the carrier.

The M&M detector is looking at the peaks, while the Gardner method is looking at the peak and zero crossing to get timing estimates.

In all cases, you'll take the error estimate, feed it into a loop filter whose output is then used to advance or retard the symbol clock for the next timing estimate. The loop filter can be thought of as either a low-pass filter, a smoothing filter or simply as weighting older estimates less than newer ones. The loop filter controls a) how fast the timing converges and b) how much noise affects the estimates.

Another option is to use an equalizer which will then remove timing delay from the signal. This can also help for very noisy signals I believe.

I have used these resources in forming this answer:
http://www.comlab.hut.fi/opetus/333/2004_2005_slides/CarrierTimingRecovery.pdf
http://mobiledevdesign.com/tutorials/radio_matched_filtering_timing/

This article is especially helpful: http://rfdesign.com/images/archive/0901Litwin32.pdf
Here is a link to Gardner's original paper from 1986: http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1096561&tag=1

Note In general this area is called "clock and data recovery" or "CDR". Specifically we're dealing with "symbol synchronization". These terms should aid you in your search.