Can I have some feedback on this research paper? I have submitted this paper "Effect of Transducer Positioning in Active Noise Control" to arXiv.

I would like to get some feedback on the methodology of the same. Inputs of any sort are welcome.
Can you please take a look at it and give me some feedback/suggestions/critics?

Comments provided here are in two broad categories: Presentation and Subject matter. The "Presentation" section is the easiest to address. There are some things that could be rephrased in terms of language use but these might be just personal preferences. The "Subject matter" section includes comments in methodology which could take more time to address.

Overall, my perception was that the paper is an excellent start to more complex work but certain key points for this type of work were not clear enough. There is always "more work" that can be done but this particular paper has room for improvement before hitting those limits. I would give it another go before submitting to a journal.

Presentation

"Passive noise control uses bulky sound absorbent materials to mas unwanted sound, whose performance...", The control's or the material's?
"In the existing ANC literature, which has a history of 7 decades, the acoustic...", do we really need to load that phrase with yet another point on the way to saying something about the existing literature? It breaks the flow and I am not sure how much it adds to the point.
"This approximation, used for ease of computation, limits researchers....". Same comment applies here. Just say what you want to say in one phrase with two strong points rather than having one main point and an aside.
"...treatment of the research work is available at [some URL]...". Why not reference this properly and include one more citation? Or add a footnote if you think it would be better to have the URL immediately visible for some reason. As a reader, I am not even going to "parse" the URL unless I want to visit it and I have not yet read the message that this paper is trying to convey. Why send the reader away already?
"...using single or multiple microphones and passes [IT] on to the...","...THE TOOLBOX provided by Habets[20] HAVE been..." I am not going to flag each and every one of these tiny little things, please proof read the document. A typo can be seen a bit like "Oh come on man! I have worked my bottom off in this paper and you are picking on a typo!?!?!". But the thing is, the typo is brakng the flw. You are nt the redr of the papr, someone else is and they have no idea what you are trying to tell them. They "form" the idea by reading these lines, word-by-word. Typos get in the way. I stop, re-read the phrase, start again...It's more effort.
"ANC PRELIMINARIES" goes on for 1.5 paragraphs BEFORE revealing that what is described is actually depicted in Fig 2. For 1.5 paragraphs, I have been building this sketch up in my head reading the description. More effort on my part (the reader...not me personally).
"...Room Acoustics Models(RAM) are often used to recreate and reproduce room acoustics without any cross-model mismatch..." What is "cross-model mismatch"?
"...beamforming, acoustic auralization, speech.." Auralization is already acoustic.
"A RIR model channel [Fig 3]...". The figure's caption could have RIR captialised "Room Impulse Response...".
Are you sure you need II and III as separate sections? I don't know the audience of the paper, but maybe these two sections are at a kind of...very introductory level (?).
You start IV with "Figure 4 shows...", that's great. Why not do the same with the rest of the figure references?
Caption of Table 1 reads "Simulation settings used for an arbitrary case". We will come back to the "arbitrary case" later on in section 2 but as far as presentation of the work is concerned "ARBITRARY CASE" reads oddly. Perhaps this is the "TEST CASE" or "SIMULATION PARAMETERS"?
Add a grid to all of your graphs. It helps with reading the graph a lot.
"The source signal considered is sinusoidal with a wideband component", "We used the same method to generate the source signal as used here[22]...". Everything that is key to the results of this paper must be sufficiently mentioned in THIS paper. I (the reader...) am not going to go to your github repository to check your code. I may have a system of my own and I am just looking to improve or test it. I would really appreciate a simple formula there showing to me how was the "sinusoidal with a wideband component" created.

Subject Matter

Section IV:

"The room size was chosen arbitrarily to a size of...". Why arbitrarily? What is wrong with "The average auditorium size is...therefore we focused on rooms of that size", or, "In our research we chose to focus on rooms of size BLAH because that's the standard listening room size". Maybe you can refer to some sort of authority that is sending standards for room sizes. Maybe the equipment you develop is ideal for rooms of BLAH size? In any case, for me "...chosen arbitrarily..." is not perceived positively. The point is why did you choose that room size?
"The z level of [THE] microphone, noise source and microphone were fixed to a height of 1.53..." WHY?. First of all, 1.53 (what? Probably meters) and then WHY?. Maybe that's the level an average person's ear ends up when they are sitting? Otherwise...WHY 1.53 meters? Is that the average size person? Is there any other reason?
"...for ease of computation....". No.
"The x coordinate[S] [...] were chosen randomly...". Why randomly? What is the reason? I don't know if it is clear what I am saying here. I am not trying to put you on the spot. What I am reading is "We took some tools, we threw some dice, we got these numbers"...Which, I don't think is what you are trying to say here. Pin these parameters down because of a reason. Then the paper gets some perspective with reality.
"The antinoise source positions were evenly spaced at a distance of 11cm". WHY?.
"Table 1 shows an arbitrary case..."... Shows one of the test cases considered, see previous comment... There is too much "arbitrary", "random".
*"The sampling frequency (Fs) is fixed to 2kHz"....WHY?. For example, you could say something like, "Given the current performance of DSP / speakers, realistic ANC is only possible for frequencies between BLAH and BLAH Hz. We targeted this band and therefore worked with Fs=2kHz.".
"...according to Nyquist theorem...". No.
*"The adaptive filter [ORDER] was fixed [AT] 350 [COEFFICIENTS] at Fs=2kHz which is an important parameter determining ANC performance."... WHY?. Why is 350 at 2kHz important for performance? Are you talking about computational performance or whether the ANC does its job?
"A step size of .." It is not entirely clear that you "sweeped" the locations of the mics and the speakers.

THE POINT ABOUT FIR FILTERS: In the introduction it is stated "In the existing ANC literature, [...], the acoustic channels are mostly modelled using [FIR] filters. This approximation, [...], LIMITS researchers in studying the effects of physical aspects like transducer positioning, orientation, room reverberation, wall reflections, room structure etc. [...] RIR is analogous to AN IMPULSE RESPONSE in LTI systems. The IR is considered to characterise the behaviour of an LTI system completely.....". Given that you are using FIR filters everywhere in the paper...Shouldn't this part be rephrased / removed completely?

The 350 tap LMS is FIR
RIRs are FIR.
Equations (1) and (2) are FIR

THE NATURAL ROOM FILTER. A room is already a filter. When you pump sound in the room, as soon as you have reflections, you get interference. Interference puts a grid of "nodal points" in the room. There are certain locations in the room that if you put a microphone there, naturally, you will be recording a very low signal. This is the reason behind positioning speakers at certain distances from the walls so that you have the waves projected at the back of the speaker work to your advantage. I am looking at your Fig 9 and I am thinking that I can almost see the grid in the Y direction. The point is: Given a room, a mic and a speaker, the mic might be recording reduced sound levels because of its natural positioning. WITHOUT ANC. The room itself is filtering. Fig 9's caption says "...It can [BE] seen that the good anti-noise source positions..."...perhaps these are not "good", but "natural" (?). Because of the physics...That's where you expect them to be. Not because of the ANC. So, what you are saying is, I am taking into account the natural acoustics of the room to correct for what the signal sounds like at the point of the mic knowing the room it propagates and this helps in ANC (?)...
Fig 11 shows that the ANC is working. This is an ANC, clearly targeting the band $(0 \ldots \approx 1000Hz]$. What is its effect on useful signal that occupies the same region?
SIMULATION METHOD / NOVELTY:. Figure 3 implies that you are using the ray tracing model of sound propagation, which, for a rectangular empty enclosed space is fine. But we are not in the 70s. So, in the presence of this fantastic (and FOSS) piece of software, maybe you could use impulse responses from realistic (and large) auditoriums with specific constraints (?). In i-Simpa you can control everything...And it's free...Which means that you can download it now and start playing.

Hope this helps.


How to remove wind noise from video file I've some video footage where wind noise is drowning out background noise I want to be clear.

Is it possible to do some DSP to remove/attenuate this wind noise in order to hear the engine noise?

I am not well versed in DSP, but my initial assumption would be;

Extract audio from video file
Calculate FFT of audio file so as to analyse the frequencies at which wind is present
Design and implement a filter which removes unwanted noise from audio file
Modify original movie file by replacing the existing audio with the "cleaned" audio

Any help appreciated.

Extracting Audio from Video

There is a simple cross-platform tool for manipulating / accessing video & audio files called ffmpeg. I've linked you the documentation available for the whole toolkit, just so you can glance around and get a feel for what sort of wonders can be performed with it.

But, all you need to do is extract audio from your video file. You can do this with ffmpeg and a pointed explanation of how to do exactly this, following this link. It's more of a tutorial / user's guide and less overwhelming than the documentation I linked above. You have two choices...

Extract the audio straight from the terminal:

$ ffmpeg -i C:/path/to/file/your-video-file.mp4 -vn -acodec copy your-audio-file.aac

Remember to change the file extensions depending on what video file extension you're working with.

Extract the audio from the terminal through Python:

import subprocess as subproc
command = "ffmpeg -i C:/path/to/file/your-video-file.mp4 -vn -acodec copy C:/path/to/your-audio-file.aac"
subproc.call(command, shell=True)

Removing Wind Noise By Analyzing FFT Data

Built-in tools for audio noise reduction / removal, like the Noise Removal Tool in Audacity, are all various implementation of the same concept: spectral subtraction. Spectral subtraction is so widely used because... it's relatively simple. All it takes is to compute the FFT of the audio signal in question, split the FFT into smaller "frames", identify the power level of the steady background noise within each frame, and then filter out that energy from the entire frame. Coming up with what filter, what frame size, and what energy level works best for your case is part of what makes each spectral subtraction technique different. So, there's no one, static solution for all audio files that I could show you - just analyze the magnitudes of your audio's FFT frame-by-frame and figure out what works best for your audio. In the end, you're left with whatever additional energy is present in the signal, which is usually the part of the interesting part of the audio.

If you don't know how to compute & view the FFT of a signal in Python, consult the Examples section of the numpy documentation on the numpy.fft.fft function.

Patching Together Reworked Audio with Original Video

Again, ``ffmpeg` is your friend here. This is as simple as a one line command:

ffmpeg -i your-video-file.mp4 -i your-audio-file.aac -c:v copy -c:a copy output.mp4

I hope that this approach helps give you the DSP experience you're looking for, while also solving the noise reduction problem you have with your current video. Good luck!


Obtaining impulse response by means of cross-correlation between input and output Having a system with an input (sound emitter) that is defined as a unit-variance white noise sequence, and I get the output (sound receiver) so I can measure cross-correlation functions between them, how could I employ the correlation function to obtain the impulse response h[n], which defines the media?

we can calculate it!

define $y=h*t$, $r_{xy} -$
is the cross correlation

$r_{xy}(\tau)=\int_{-\infty}^ \infty x(t)y(\tau+t)dt = (y*x^-)(\tau) = (x^-*h*x)(\tau) =(r_{xx}*h)(\tau)\underbrace{=}_{r_{xx}=\delta}h(\tau)$

so actually the cross correlation is the impulse response


What are some good images to use when comparing image compression algorithms? There was a question about images that are hard to compress. Are there images that are better suited to testing compression quality than others?

Maybe an image of diagonal lines of specific colors at a specific angle would be more likely to show artifacts of algorithms / libraries than others? High contrast colors? High contrast black and white? Checkerboards? Certain kinds of photos?  Maybe there is a set of standard images to use to compare? Is it mostly subjective or are there objective criteria for comparing?

Are there images that are better suited to testing compression quality than others?

You are exactly right. A good compression algorithm is one that performs well on average considering all those types of images. In reality, a database of different types of image is used to evaluate a particular compression techniques.

  Maybe an image of diagonal lines of specific colors at a specific angle would be more likely to show artifacts of algorithms / libraries than others?

Good test image makes sense when it comes to a specific feature of image we are interested to look at. For instance, assume we have an image compression based on discrete wavelet transform. Now we are interested to see how much gain can be achieved in terms of compression ratio if we replace the ordinary wavelet we use with a directional wavelet. In such case, an image full of straight lines with various directions can reveal the differences. An example test image of this type is Barbara. Or we may also use an image with a lot of smooth curves to see how the directional wavelet can decompose the curves. Lena image has more curves in this case.

  High contrast colors? High contrast black and white?

Let's say we want to study the impact of smoothness of the wavelets we select on the compression ratio. So we need two groups of test image: highly edgy images with a lot of discontinuities where the non-smooth wavelet will perform better, and low-contrast images where the smooth wavelets are more appropriate.

Or assume the cross-channel compression capability (e.g. dependencies between chroma and luma) is the subject of test. So we need to test images that have (or don't have) significant dependency across these channels.

  Checkerboards? Certain kinds of photos?

Sometimes block artifacts need to be tested considering different block sizes or at different scales. Or sometimes the compression algorithm is designed for a specific purpose such as remote sensing or satellite images. So we need to exclude e.g. portraits and focus on aerial test images. Or similarly for medical imaging...

  Maybe there is a set of standard images to use to compare?

There are some known test images such as Lena, Barbara, Cameraman, etc. But they are mostly for academic papers, since their metrics (like PSNR at different BPPs) are quite widely known and people can evaluate the proposed techniques by looking at the reported metrics. Something like: "Lena (512) with 5-level cdf 9/7 with SPIHT at 0.5 BPP and arithmetic coding can achieve around 37.2 dB PSNR".

However, in practice (for design of a codec) a huge database is used. You can see for example something like here, here, or here.

  Is it mostly subjective or are there objective criteria for comparing?

So based on the above discussion, we can see both subjective and objective tests apply, depending on what we are interested in. The compression performance of the image coded in general, or are we just want to test a specific idea or a particular part of the whole compression pipeline.


Improving Velocity Estimation Using Multiple Sensors in a Dynamic System I have a sensor reduction model which gives me a velocity estimate of a suspension system(velocity 1) .

This suspension system estimate velocity is used to calculate another velocity(velocity 2) via a transfer function/plant model.

Can I use velocity 2 to improve my velocity estimate (velocity 1) through Kalman filtering or through some feedback system.??

You can do that if you have new data to merge in.
For instance if you have a sensor for Velocity 1 and a sensor for Velocity 2 and in addition a connection between Velocity 1 and Velocity 2 you can use Kalman.

Yet if you have only one sensor for Velocity 1 which with the transfer system you estimate the Velocity 2 there is no additional info to fuse in.

What might be done is that if you have information about the model of Velocity 2 in time you can use that (Though, the right thing to do would be reflect this model into Velocity 1 by the connection of the two, apply Kalman on Velocity 1 and from there infer Velocity 2).


Side lobes not visible in FFT of rect function I was making a plot to show side lobes for different spectral windows and I noticed if I take a simple fft of a rectangle window, I don't see any side lobes. I think I'm just making a simple mistake.

Here is the Python code I ran:

from matplotlib import pyplot as plt
import numpy as np
fAxis = np.linspace(-1/2,1/2,100)
rect = np.ones( 100 )
specWindow = abs(np.fft.fftshift(np.fft.fft(rect)))
plt.plot(specWindow)
plt.show()

Could anyone explain how to change my code slightly to see the side lobes I'd expect. I expect a discrete FT of a rectangle function to give me the Dirichlet kernel.

You have to imagine your signal as being periodically extended forever in both directions. Since you square window is constant over the whole interval when it is periodically extended it becomes flat everywhere. In other words, you don't have a square window. You just have a flat line. That's why only the zero frequency band on the fft has any amplitude. Try doing it again with some zeros on either side of your square window.


Parallel Processing Gigapixel Images Parallel programming is often used to accelerate image processing methods (e.g. Eklund et al., 2013); however, there doesn't appear to be the same level of discussion about using parallel processing to handle gigapixel images (or images that are too large to load on one machine). Is there a reference that discusses parallel processing techniques for large images?

Follow-up Question: To me, it seems there is only one way to process gigapixel images: break them up into chunks and process each chunk on a different node. The way you split the image could differ, but the basic procedure stays the same. Are there other ways parallel processing can be used to handle gigapixel images?

Parallel processing of large images or large quantities of images is no different than parallel processing for other data types. You will come across the same strengths and limitations.

First of all, there are two major paradigms for parallel problems: Those that are obviously parallelisable and those that are not (or inherently serial).

In an obvious parallel algorithm, the result of a computation can be achieved by the repeated application of one or more functions (the function here as an element of computation) to a dataset. But, in a non-obvious parallel algorithm, the result of the computation cannot be "untangled" and has to be carried out as a single function. There are of course intermediate cases with a task being composed of some parts that can be parallelised in a straightforward way and some tasks carried out serially (unavoidably).

If you are looking for an object in an image through cross correlation, you can apply an image pyramid on it and then feed the different scales to different algorithms that slide the template over the image and try to find interesting points.

If you are trying to derive the "map" of a metric, again, this can be obviously parallelised. For instance, generate an image where each pixel is the standard deviation of the 8 surrounding pixels in the original image.

If you are trying to extract SIFT descriptors, again, this can be parallelised, to a certain extent.

If you are trying to stitch images together, what you cannot parallelise easily is the optimisation process of trying to match the locations of the keypoints. This is because, you are trying to minimise an error that is distributed over many different combinations of keypoints. Therefore, inevitably, even if you tried to do it path-by-patch you would still have to tune the error between patches to make sure that the algorithm has converged to a single optimal solution.

Given what is mentioned in the comments of the question, it might be possible to still converge to a solution as far as the keypoints are concerned. But, re-projecting already large images and then putting them together into one ultra huge image, might still be challenging, not because of lack of parallelism but because of lack of memory.

For more information, please see this link and this link

Hope this helps.

EDIT:

  In the simplest form, I'm attempting to register two gigapixel images A and B

Gigapixel images will require a huge amount of memory if processed in a "conventional" way. This is more or less a fact. The question now is what do we do about it?

I say "more or less" because there are techniques that rely both on software and hardware that could allow in memory processing of such images in the "conventional" way. These are in-memory computing and Shared (or Distributed) Memory. On the former, I forget the exact term right now, but there is a technique which links the actual physical memory of a number of machines together through network links at a very low level. So, suddenly, an algorithm has huge amounts of memory available as if it was "local".

But short of these little beasts the other thing that one can do is process the images in an "unconventional" way. And in this case it really depends on what you want to achieve.

The trick is similar to Sparse Matrices.

Sparse matrices present an interface to the user that is no different than a conventional matrix. You still say I[1200,1500] but internally, the getitem operation of I's data type, that is implied by the use of indexing, does a bit of "magic" with the indices and returns a value. If a matrix is 100x100 and has at most 45 non zero double precision floating point values in it, there is no reason to waste the rest of the space. So, internally, when I[20,32] = 3.1415928 you might store something like 20, 32, [data value] and then later, when 2 * I[20,32], you look up the indices and if they exist in your lookup table you return the value but if not, then you return a 0.0.

In the case of huge images, we are on the other end of the spectrum but still the solution is kind of similar. You create something like a class MappedImage. Instead of loading the image, you have a constructor like MappedImage(pathToDirectory). The contstructor loads the image tiles and does the cross registration but does NOT! actually cross register the images. In other words, all you are after at that point is the transform that takes you from "image space" to "cross registered space".

Once this MappedImage is constructed, you still access it via I[789,145890987] (I here being an object of MappedImage) but the getitem of MappedImage now takes the desired indices, passes them through the transformation matrix, gets the actual locations, "guesses" which image patch contains the pixel values and reads the actual pixel value from that image only.

In this way, you don't have to be keeping all images in memory all the time and you can still use standard memory mapped I/O through the facilities of your language of choice.

The paper provided is primarily focused on doing a high quality registration not necessarily accessing the cross registered image. In a way it does highlight the same problems highlighted earlier and hints at the "conventional" solution which is "work on many small images".

Now, obviously, the next question is "What do you want to do with this 'image'"?

Because, if you want to visualise it, then you can create a MappedImage_Tiler that accesses some I via the interface mentioned earlier and cuts this huge image that you will never be able to see in its entirety into smaller pieces. These can then be sent to an off the shelf tile image viewer and allow you to browse the image like (for example) Google Earth.

But, if you want to work on it and produce some directional derivative or whatever else, then, yes the I[34,2875623897652398756] interface will still work but you will have to lower your expectations on speed. Especially if the goal is to work on "laptop type" hardware as the constraint was in the paper cited.

  I can then map those transformation to higher-resolution levels; however, there is error in this type of mapped registration. I'm hoping to use parallel techniques to register A and B directly at the highest resolution.

This is what I had in mind too: Use the low resolution image derived transform and "project" it to the high resolution images (possibly just through a magnification component (?)).

Hope this helps.


How can an FFT plot not be symmetrical? From discussions all over the internet, I understand that in a 512 sample point fft, only the first 256 contain relevant information provided that my signal is real. The other half contains complex conjugates of the first half. I base my question on this discussion at stack exchange.

I'm trying to plot the fft of a wav file. I used 2 files - one wav file containing some music and another that was generated by scrip to play at 440Hz.
The fft plot of the 2nd (440hz signal) showed 2 peaks - one at 440hz and another at the other end of the plot. This was expected. However, plot of the music wav file was asymmetric. How can that be?

I used the file onclassical_demo_elysium_anonymous-elysium_the-young-false-man_small-version_live-and_restored found here. Scipy did not parse it properly first and I had to strip the meta tags using audacity before I could get it working.

I used print  and I got this for the 440hz tone :

[ 203395.00000000  +0.j          203395.08226712 -57.78564066j
203395.32901382-115.57039787j ...,  203395.74007235+173.35334037j
203395.32901385+115.57039784j  203395.08226713 +57.78564063j]

For the music file, I got this :

[[ 2.+0.j -4.+0.j]
[-1.+0.j -7.+0.j]
[ 0.+0.j -2.+0.j]
...,
[ 0.+0.j -6.+0.j]
[ 0.+0.j -2.+0.j]
[ 1.+0.j -1.+0.j]]

I don't think any of them are symmetrical! Or else the first and last values should coincide at least in the real part right?

And just for the sake of it, the 440hz signal plot :

I tried plotting abs(fft) - still asymmetrical. Also, why is the plot having both blue and green colors - real and complex parts?

Using the wav version of the file you linked to I performed an fft in matlab and got:

First 3 terms:

 left channel       right channel
-1.1594 + 0.0000i  -1.8450 + 0.0000i
-0.4186 - 1.0013i  -0.2851 - 1.7655i
 0.5507 - 0.7589i   1.0478 - 0.4655i

Last 3 terms:

 left channel       right channel
 0.7388 - 0.3032i   0.4740 - 0.6654i
 0.5507 + 0.7589i   1.0478 + 0.4655i
-0.4186 + 1.0013i  -0.2851 + 1.7655i

Looks conjugate symmetrical

I wonder if you are processing the left and right channels at the same time and it is thinking it is a 2d signal or a complex valued signal.

The only way you get an asymmetric fft result is if the signal is complex.


Confusion about how to calculate signal to noise ratio(SNR)? I searched a lot but I didn't find a good answer. According to wikipedia SNR page, SNR can be calculated by both ratio of Power of clean signal to its noise Power either ratio of variance of clean signal to variance of its noise, but actually power and variance are two different things so I am confused how we can use both of them equally. I should mention that I speak about two finite signal, first the original one and second the noisy one and want to calculate the SNR.
Thank you

The SNR is the ratio of the powers. Note that for zero-mean signals, the variance equals the power, so the SNR is the ratio of the variances of signal and noise, but only if both are zero-mean.


U-nets : how exaclty is upsampling performed In U-nets I would like a more straight-forward/detailed explanation in how the upsampling part ("right part of the U") is performed.

I read that it can be done by "transposed convolution layers" aka. "deconv layers". I would like a clear example (possibly a bit simplified) on how that is performed?

e.g. initialization, why/how? Detailed sizes/shapes input output etc?

e.g. from here https://medium.com/@keremturgutlu/semantic-segmentation-u-net-part-1-d8d6f6005066 :
"In transposed convolutions we have weights that we learn through back-propagation. In papers I’ve come across all of these upsampling methods for various cases and also in practice you may change your architecture and try all of them to see which works best for your own problem. I personally prefer transposed convolutions since we have more control over it but you may go for bilinear interpolation or nearest neighbor for simplicity as well."

I don't understand the stated "ou may go for bilinear interp etc."

why would we choose a fixed interp method over learned transposed conv2D filter?
How can that make sense?

The original U-Net paper used nearest neighbor interpolation as far as I know. This is also the default upsampling method in TensorFlow. My own anecdotal advice is to not use transposed convolutions in U-Net. It will only make your CNN slower and won't really increase your F1 score (or whatever metric you're using). I would recommend nearest neighbor interpolation.
Transposed convolution is basically regular convolution with zeros inserted between each input. In U-net after each upsampling, there will be a regular convolution layer. So you don't lose anything if you choose fixed weights like nearest neighbor.
Upsampling outside the neural network domain applies two steps: 1. upsample, 2. lowpass filter. The lowpass filter is here the learnable convolution. In CNN there only seems to be an issue with aliasing (see Making Convolutional Networks Shift-Invariant Again). But this happens also with max pooling / downsampling.


Determining effect of powers of DFT matrix conjugate on input without multiplying Suppose I have a DFT matrix $F$ of dimensions $N \times N$ where $N$ is a power of $2$. I convert it into a unitary matrix $\displaystyle U = \frac{F}{\sqrt{N}}$ and compute $S = U^* U^* $.

$S$ seems to have a pattern. $S$ has at most one $1$ per row. Applying $S$ on my $N$-sample signal $x$ essentially provides me either a permutation of $x$(if $N$ is an even power of $2$) or a subset of $x$(if $N$ is odd power of $2$). By subset, the length of $x$ remains the same, but some of the elements get truncated to $0$

Translated to MATLAB code for $N=8$,

N=8;
F = dftmtx(N);
U = F/sqrt(N);
S=ctranspose(U)*ctranspose(U);

$S$ above looks like the following:

My problem is the following: Given $N$, can I determine how to pick elements of $x$ without having to resort to multiplying by $S$ ? (For the curious, let's say the matrix multiplication is expensive for my situation).

In both cases (odd and even powers of 2) you should get as a result simply a time inversion of your original signal. Because of the implicit periodicity you get from the original vector

$$[x_1,x_2,x_3,\ldots,x_N]$$

the following time-inverted (and periodically continued) vector

$$[x_1,x_N,x_{N-1},\ldots,x_2]$$

Here's an 'odd power' example:

N=8;F=dftmtx(N);U=F/sqrt(N);
x=[1 2 3 4 5 6 7 8]';
y=U'*U'*x;
y =

1.00000 - 0.00000i
8.00000 + 0.00000i
7.00000 + 0.00000i
6.00000 - 0.00000i
5.00000 + 0.00000i
4.00000 + 0.00000i
3.00000 + 0.00000i
2.00000 - 0.00000i


ARMA & MA methods: how do you know the error terms? Reading the ARMA model for the first time, and I'm confused.

Let's say I have a time series

x = [1, 2.1, 2.9, 3, 4.1]

According to the ARMA model, $X_t$ is a linear combination of previous values and errors, something like

$$\sum_i \phi_i X_{t-i} + \sum_i \theta_i \epsilon_{t-i}$$

But, what are $X_i$ and $\epsilon_i$??

is $X_i$ the actual value of the series at $t=i$? Eg, in my example $X_2 = 2.1 ?? (with 1-based indexing)
if so, what are the error terms?

The same question for the simple moving-average model, where all the $X_i$ values are remain unused.

In the context of discrete-time statistical signal processing, an ARMA-(p,q) random process is defined as (assuming zero initial conditions)

$$ \sum_{k=0}^p a_k ~ x[n-k]  = \sum_{k=0}^q b_k ~ v[n-k] $$

or equivalently

$$ x[n] = - \sum_{i=1}^p a_i ~ x[n-i]  + \sum_{i=0}^q b_i ~ v[n-i] $$

where $v[n]$ is a white-noise (WSS) random process with variance $\sigma_v^2$ and $x[n]$ represents the resulting ARMA process.

In this context your coefficients are related as $$\phi_i = - a_i $$ and $$\theta_i = b_i$$ The process reggresses over its past values $x[n-i]$ ($X_{t-i}$ acc) and is also a moving average of the input noise $v[n-i]$ ($\epsilon_{t-i}$ acc).


Is there an audio renderer? I wanted to know if there is a sound render like the renderers of 3d modeling software (vray ,  cycles etc).

I would like to experiment with sound in a controlled environment

if it does not exist which is what makes it difficult to create an audio renderer?

"Sound renderers" do exist, and some are very good: have a look at the research of Nikunj Raghuvanshi. For an introduction to their rendering method:

Nikunj Raghuvanshi, Rahul Narain and Ming C. Lin, Efficient and Accurate Sound Propagation Using Adaptive Rectangular Decomposition, IEEE Transactions on Visualization and Computer Graphics(TVCG), 15(5), 2009.

Abstract:

  Accurate sound rendering can add significant realism to complement
  visual display in interactive applications, as well as facilitate
  acoustic predictions for many engineering applications, like accurate
  acoustic analysis for architectural design. Numerical simulation
  can provide this realism most naturally by modeling the underlying
  physics of wave propagation. However, wave simulation has
  traditionally posed a tough computational challenge. In this paper, we
  present a technique which relies on an adaptive rectangular
  decomposition of 3D scenes to enable efficient and accurate simulation
  of sound propagation in complex virtual environments. It exploits the
  known analytical solution of the Wave Equation in rectangular domains,
  and utilizes efficient implementation of Discrete Cosine Transform on
  the GPU to achieve at least a hundred-fold performance gain compared
  to a standard Finite Difference Time Domain (FDTD) implementation with
  comparable accuracy, while also being an order of magnitude more
  memory-efficient. Consequently, we are able to perform accurate
  numerical acoustic simulation on large, complex scenes in the
  kilohertz range. To the best of our knowledge, it was not previously
  possible to perform such simulations on a desktop computer. Our work
  thus enables acoustic analysis on large scenes and auditory display
  for complex virtual environments on commodity hardware.

The paper has its own website. The manuscript version of the paper there is missing some prose, so preferably get the paper from the first link.

Siemens, for example, advertises acoustics simulation as part of their computer aided engineering solutions.


Audio signal processing book suggestions I want to start learning about audio signal processing but I cannot find a good starting point. Please suggest books or lecture series.

If you can get it:

Zölzer, Udo (1997). Digital Audio Signal Processing. John Wiley and Sons. ISBN 0-471-97226-6

or

Orfanidis

might have to buy them used.

is your interest music processing and synthesis?

there are some nice online books.  Look for anything Julius Smith.  But there are others, but i can't remember any names to search for.

Remembered one:

Davide Rocchesso

here's another link

another old one:

Ken Steiglitz


Parseval's theorem - Average Power I want to calculate averaged power of a time-domain signal by means of its spectrum. I guess Parseval is the right tool.

So I sample a sinus of 100 Hz 10000x within one second.

Unfortunately the sum of the squared samples euqals not the sum of the FFT amplitudes (weighted by the number of FFT bins). Where is the mistake?

# Some python code

import numpy
import matplotlib.pyplot as plt

# Create Time Domain Signal for 1 sec
fs = n = 10000    # Samplingfrequency
ti = numpy.linspace(0,1,num=fs)
sx = 1*numpy.sin(2*numpy.pi*100*ti)

# Calculate spectrum via FFT and account for scaling n/2
# taking the real fft (rfft) only the positive frequencies are calculated
fx = numpy.fft.rfft(sx)/(n/2)
no_of_points = fx.shape[0]

# Calculate RMS for time domains signal + spectrum
parseval_sx = numpy.sum(sx**2)
parseval_fx = numpy.sum(numpy.abs(fx)**2)/no_of_points

print parseval_sx, " equals not ", parseval_fx

Output:

4999.5  equals not  0.000199940012002

It's possible to calculate average power (proportional to RMS amplitude) with the rfft, but it's more complicated because the DC and Nyquist bins are not doubled in the full fft (and the Nyquist bin doesn't exist if n is odd).

NumPy's rfft is preferred over SciPy's, because it outputs complex numbers, which are compatible with N-dimensional arrays, rather than SciPy's "packed" format.

def rms_rfft(spectrum, n=None):
    """
    Use Parseval's theorem to find the RMS value of an even-length signal
    from its rfft, without wasting time doing an inverse real FFT.

    spectrum is produced as spectrum = numpy.fft.rfft(signal)

    For a signal x with an even number of samples, these should produce the
    same result, to within numerical accuracy:

    rms_flat(x) ~= rms_rfft(rfft(x))

    If len(x) is odd, n must be included, or the result will only be
    approximate, due to the ambiguity of rfft for odd lengths.
    """
    if n is None:
        n = (len(spectrum) - 1) * 2
    sq = real(spectrum * conj(spectrum))
    if n % 2:  # odd
        mean = (sq[0] + 2*sum(sq[1:])           )/n
    else:
        mean = (sq[0] + 2*sum(sq[1:-1]) + sq[-1])/n
    root = sqrt(mean)
    return root/sqrt(n)

More details here: https://gist.github.com/endolith/1257010


Epoched data to a continuous signal if we create a continuous signal from epoched data (say, EEG data with 100 epoches), what kind of problems it may cause( i am only aware of edge artifacts ) and how i can avoid them in analysis (say, power analysis)? Thanks in advance,

You may have problems related to discontinuity and this can give small glitches in your signal, if you hear your signal you will notice pops/clicks artifacts at every 100 epoches edge, one way to try avoid this problem is apply one window function in a overlapped segments between your Epoches.

When you apply a window like a Hann window this reduces the amplitude of the components of the outlet ends of the window and its suppress the edge artifacts when reconstructing your signal to be continuos, for it you need use technics like overlap and add to works like a crossfading between every overlapped segments.

OK you asked for example, i can show you the basic concept, imagine two sinusoids at 100hz, if you try join you will have problems, take a look:

In this plot is clear that the end of the sine wave does not match with the beginning of another, and this cause the problem described above.

One way to do try avoid it is apply a crossfade, here an example:

Fs=22050;
f=100
winsize=2048;

%creating signal
signal1= 0.9*sin(2*pi*f/Fs*(1:winsize));
signal2= 0.9*sin(2*pi*f/Fs*(1:winsize));

%use plot(in) to see the first image plot in this answer
in=[signal1 signal2];

%create a hann window
window1=hann(winsize)';
window2=hann(winsize)';

%i'm using just the last half values of the hann window to skip fade in
window1(1:winsize/2)=ones;

%i'm using just the first half values of the hann window to skip fade out
window2(winsize-(winsize/2):winsize)=ones;

%apply window in signal1 and signal2
teste1=(signal1.*window1);
teste2=(signal2.*window2);

%apply 50% crossfade between signal1 and signal2
OverlapCrossfade=teste1(length(teste1)/2:length(teste1)-1) +  teste2(1:length(teste2)/2);

%rebuild the signal after crossfade
join = [teste1(1:length(teste1)/2) OverlapCrossfade teste2(length(teste2)/2:length(teste2)-1)];

Now if you do plot(join) you can see it:

OW this is better i cant see glitches (but we can see amplitude variation), now all that you need to do is apply some type of envelope adjust.


What is spectral entropy? Can anyone explain what spectral entropy is?

Does noise with a restricted bandwidth have the same spectral entropy as white noise?

Spectral Entropy describes the complexity of a system. It is defined as follows:

Calculate the spectrum $X(\omega_i)$ of your signal.
Calculate the Power Spectral Density of your signal via squaring its amplitude and normalizing by the number of bins.

$$P(\omega_i)=\dfrac{1}{N}\left|X(\omega_i) \right|^2 $$

Normalize the calculated PSD so that it can be viewed as a Probability Density Function (integral is equal to 1).

$$p_i=\dfrac{P(\omega_i)}{\sum_iP(\omega_i)} $$

The Power Spectral entropy can be now calculated using a standard formula for an entropy calculation.

$$PSE = -\sum_{i=1}^np_i\ln p_i $$

In case of boosting of your noise signal, without performing any other processing, the Entropy will change. I guess there is no other way around that.


Remove Shadows from Contour of an Object (Grapes) I want to extract only the grapes from the images. Unfortunately, sometimes I don't get the grapes exactly.

I have a code that doesn't work well in all cases. Sometimes, it does not differentiate between the shadows and the grapes.

The example input image:

The result I'm getting:

And this is my code:

RGB = imread('DSC02807.JPG');
GRAY = rgb2gray(RGB);

threshold = graythresh(GRAY);
originalImage = im2bw(GRAY, threshold);

originalImage = bwareaopen(originalImage,250);

SE = strel('disk',10);
IM2 = imclose(originalImage,SE);

originalImage = IM2;

imshow(originalImage);

This is what your code produces (using colormap(hot) with white boundaries that are obtained from IM2):

I think the shadows are nicely treated, why are you thinking that your code doesn't notice between shadow to grapes?

This is the code used for producing this image

bw=IM2<1;
[B,L] = bwboundaries(bw,'noholes'); % fill any holes, so that regionprops can be used to estimate
imagesc(GRAY); hold on
for k = 1:length(B)
    boundary = B{k};
    plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 1); hold on
end
colormap(hot(256))


Solving a second order difference equation I need to solve the following 2nd order difference equation

$$3y[n]-y[n-1]+y[n-2]={\bigg(\frac{1}{2}\bigg)}^{n} u[n]$$

where $u[n]$ is Heaviside's step function. I started by trying to find the homogeneous solution:

$$3-z^{-1}+z^{-2}=3-\frac{1}{z}+\frac{1}{z^2}=3z^2-z+1$$

Homogeneous equation:

$$3z^2-z+1=0\leftrightarrow z=\frac{1\pm\sqrt{1-4\times 3\times 1}}{2\times 3}\leftrightarrow z_1=\frac{1+\sqrt{11}i}{6};z_2=\frac{1-\sqrt{11}i}{6}$$

Since I have complex numbers I don't know how to write the homogeneous solution with the constants $C_1$ and $C_2$. I also don't know how to find the particular solution.

You can indeed proceed with the complex roots $z_1$ and $z_2$ as well as with the real ones.

In order for a particular solution you may assume a solution of type:
$$y_p[n] = K ( \frac{1}{2})^n u[n]$$

And plug it into the LCCDE to find the undetermined coefficient $K=1/5$.

Then assuming a homogeneous solution of the form (for causal system) $$y_h[n] = C_1 z_1^n u[n] + C_2 z_2^n u[n]$$

You have the complete solution as:
$$y[n] = y_h[n] + y_p[n] = (C_1 z_1^n + C_2 z_2^n + \frac{1}{5} (\frac{1}{2})^n) u[n] $$

In order to find $C_1$ and $C_2$ you can get $2$ equations from $y[0]$ and $y[1]$ by explicitly solving them from the LCCDE via recursion, which yields
$y[0] = 1/3$ and $y[1] = 5/18$. Then you have the following:

$$
\begin{cases}
y[0] = \frac{1}{3} &= C_1 + C_2 + \frac{1}{5} \\
y[1] = \frac{5}{18} &= C_1 z_1 + C_2 z_2 + \frac{1}{10} \\
\end{cases}
$$

Solve this system for $C_1$ and $C_2$ to complete the procedure...

Alternateley you can also consider the following.
When have a LCCDE with real coefficients, its consequence is the fact that the roots associated with the homogeneous equation will be either real or in complex-conjugate pairs, as you have found.

Furthermore you can deduce from the LCDDE with real coefficients that its solution $y[n]$ must be real valued, for a real input $x[n]$.

Hence you can deduce from these that the constants $C_1$ and $C_2$ of those two complex roots should be constrained like $C = C_1 = C_2$ or $C = C_1 = C_2^*$  then you can see that the complex roots will merge:

$$ C_1 z_1^n + C_2 z_2^n = C (z_1^n + (z_1^n)^* ) = 2C \mathcal{Re} \{z_1^n\}$$

or like

$$ C_1 z_1^n + C_2 z_2^n = C z_1^n + (C_1z_1^n)^*  = 2 \mathcal{Re} \{C z_1^n\}$$

Conversion to polar form helps like:
$$z_1 = |z_1| e^{j\phi_1} \Rightarrow 2C \mathcal{Re}\{z_1^n\} =2C |z_1|^n\cos(\phi_1 n) $$

or like
$$C = |C|e^{j\theta_C}, ~~z_1 = |z_1| e^{j\phi_1} \longrightarrow 2 \mathcal{Re}\{C z_1^n\} =2 |C| |z_1|^n\cos(\phi_1 n + \theta_C) $$


Goerzel algorithm - increasing accuracy I have an audio buffer with 1024 samples per second.
The sample rate is 44100.
I am trying to identify the frequency tones.
Using the Gorzel algorithm, I can get the frequency but it's not that accurate.
I was trying some other scope program and I got better results with the same hardware.
How can I decrease the id-bandwidth (to be sharper) so the id will be accurate for the target frequency?
(I can't use FFT because it's real time and takes some memory.)

float goertzel_mag(int16_t* data ,int SAMPLING_RATE ,double TARGET_FREQUENCY,int numSamples )
{
    int     k,i;
    float   floatnumSamples;
    float   omega,sine,cosine,coeff,q0,q1,q2,magnitude,real,imag;
    float   scalingFactor = numSamples / 2; // -2

    floatnumSamples = (float) numSamples;
    k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
    omega = (2.0 * M_PI * k) / floatnumSamples;
    sine = sin(omega);
    cosine = cos(omega);
    coeff = 2.0 * cosine;
    q0=0;
    q1=0;
    q2=0;

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

    real = (q1 - q2 * cosine) / scalingFactor;
    imag = (q2 * sine) / scalingFactor;

    //double theta = atan2 ( imag, real); //PHASE
    magnitude = sqrtf(real*real + imag*imag);
    return magnitude;
}

There are two things you can do:

Use a different algorithm, like autocorrelation.
Use the phase information! This technique is discussed in chapter 9 of DAFX. You can find matlab source here.

DAFX discussed the phase-based technique using the FFT, but it should apply to Goerzel's algorithm as well, assuming you have real data at that frequency. The key concept is that you must note how much the phase changes over time to determine the exact frequency.

You may want to experiment with an FFT of smaller size instead of Goerzel's algorithm to validate the assumptions.


High Order Filtering Which window has better performance for higher order filtering?

A. Rectangular window
B. Hamming window
C. Kaiser Window
D. Hann Window

Can anyone please answer this with explanation.

So, first of all, we need to be clear that you mean FIR filtering not IIR.

By "high-order filtering", I presume you mean an FIR with a long length:  $L \gg 1$.  That is not unusual.  FIR filters have to be pretty long in order to get anything useful done.  The sharper the selectivity of the filter, the larger $L$ must be.

Most often, the best practice for using the "windowing method" of FIR design, is to use a Kaiser Window because it's optimized to give you, for a fixed $L$, the best stopband rejection or the sharpest transition-band region (what goes in between the passband and the stopband).  There is this $\beta$ parameter that let's you adjust that tradeoff.

So, if you're doing the windowing method (say, using MATLAB), first choose and FFT length, $N$, that is much much larger than your expected FIR length.  Like your FIR length might be $L=1024$, but then set $N=65536$ or something.  Then draw out your target frequency response and inverse FFT that to see what the impulse response might look like.  Remember the latter half is really negative frequency before the iFFT and negative time after the iFFT, so you may have to use fftshift() to swap the halves.  Then apply the window to that longer FIR to make it shorter, but for the purpose of analysis, keep the zeroed samples that became zero because of the windowing.  Then FFT back to see how bad the windowing affected your frequency response.  Adjust the Kaiser $\beta$ to see what your tradeoff looks like.

Rinse and repeat.  If you can't get the frequency response you want, you may have to increase the FIR length $L$.


How to transmit voice signal and receive on radio? I want to transmit data like radio stations do.

Like radio stations broadcast music and etc. using frequency, so I am also curious to know about it that actually how it works.

I also want to make such a device that will broadcast my data or voice in my apartment and if anybody will on that frequency they should hear my broadcast data.

Is it possible or not ?
If yes then anyone please guide me.
And what happen if i ll broadcast on same frequency on which already data is broadcasting ?

I want to know about radio signals, so any guidance on it.
Thank You.

The first thing you have to know is that transmitting on the FM band is something heavily regulated. Your emitter must not exceed a specific power or range. Check the legislation of the country in which the device operates.

For this you will need an FM transmitter. You can buy them pre-made, or build them with very few parts. It might even be possible to get it to work without additional hardware and some creative hacking (some smartphones or tablets use the BCM4329 chip for example, which is capable of FM transmission).

  And what happen if i ll broadcast on same frequency on which already data is broadcasting ?

The strongest signal wins.


What exactly is transmit power in wireless communication? My background is computer science. Recently, I had to come across some concepts in wireless communication/cellular mobile communication. Despite having gone through many materials, I am still having a hard time understanding some basic concepts, especially the concept of "transmit power".

So, what exactly is transmit power? Is it:

the electrical power needed by a transmitter to transmit a signal to a receiver?
or the amount of electrical power transmitted from a transmitter to a receiver via a channel?
or is it a different kind of power?

The transmit power for a single output transmitter typically refers to the output of the transmitter power amplifier that is connected to the transmit antenna, and that power is almost entirely radiated over the air. It is power in Watts.

The power needed by the transmitter includes this and everything else to run the transmitter. The power amplifier (PA) often dominates this power requirement which is why PA efficiency is such a big deal (For example, if a PA has a 50% efficiency and is transmitting 2W, it will require 4W of DC power to run). PA's are more efficient the closer we run them into saturation, hence constant envelope modulations such as FM (which can be saturated with no further distortion) are attractive when power efficiency is of prime importance.


Matlab video processing of Contracting tissues. how to analyze movement? I'm trying to write a code The helps me in my biology work.
Concept of code is to analyze a video file of contracting cells in a tissue

Example 1

Example 2: youtube.com/watch?v=uG_WOdGw6Rk

And plot out the following:

Count of beats per min.
Strenght of Beat
Regularity of beating

And so i wrote a Matlab code that would loop through a video and compare each frame vs the one that follow it, and see if there was any changes in frames and plot these changes on a curve.

Example of My code Results

Core of Current code i wrote:

for i=2:totalframes
        compared=read(vidObj,i);
        ref=rgb2gray(compared);%% convert to gray
        level=graythresh(ref);%% calculate threshold
        compared=im2bw(compared,level);%% convert to binary
        differ=sum(sum(imabsdiff(vid,compared))); %% get sum of difference between 2 frames
        if (differ ~=0) && (any(amp==differ)==0) %%0 is = no change happened so i dont wana record that !
            amp(end+1)=differ;  % save difference to array amp wi
            time(end+1)=i/framerate; %save to time array with sec's, used another array so i can filter both later.
            vid=compared; %% save current frame as refrence to compare the next frame against.
        end
end
figure,plot(amp,time);

=====================

So thats my code, but is there a way i can improve it so i can get better results ?

because i get fealing that imabsdiff is not exactly what i should use because my video contain alot of noise and that affect my results alot, and i think all my amp data is actually faked !

Also i actually can only extract beating rate out of this, by counting peaks, but how can i improve my code to be able to get all required data out of it ??

thanks also really appreciate your help, this is a small portion of code, if u need more info please let me know.
thanks

You could try a simple non-linear denoising filter (3x3 or 5x5 median filter) on your image sequence to make your detection more robust to noise ; but I get the feeling that the moving details on these images or quite thin and that too much filtering would harm your detection. And it looks like you are already thresholding the images before comparison, so the only noise that matters here are pixels moving above and below the threshold at each frame.

To get a beating rate you could also compute the auto-correlation of the sequence. It will show a secondary maximum at the beating rate. The autocorrelation is actually an interesting representation for your task because:

A very regular beat will mean a very sharp secondary peak on the AC sequence ; while a more irregular beat will show a flattened, spread out secondary peak in the AC sequence.
The ratio of the autocorrelation secondary peak to the first peak ($r_0$) can be used as a good measure of beat strength.

But back in time domain you could also consider:

The kurtosis of the sequence as a measure of "peakedness" (and thus beat strength). Other measures of peakedness include the ratio of geometric mean to arithmetic mean.
The standard deviation of the time intervals between peaks as a measure of beat regularity. The downside is that you will have to detect peaks, which will be a bit imprecise in your case especially given the low temporal resolution of your data (and this is why I thought about doing things from the AC function first!)


DSP for computer engineers DSP is more kindred with electrical engineering than with computer science, how ever there are topics of DSP like MP3 encoding, JPEG compression and sound analysis can be understood by a computer engineer (with a previous strong math background in Z transform and Fourier series).

For a BS in computer science that only have very basic knowledge in electronic, is there a way to name that kind of specialization in DSP? I.e. someone who is looking to create DSP algorithms that can be run from a computer but not looking to implement them in a custom piece of hardware.

My question is mostly related for a computer science student who wants to apply for a master degree in an international university (no knowledge in analogic circuits, hardware description languages or microprocessors)

Would applied math, be a better approach??

DSP is a domain mixed between EE/CE and Math/CS. Hardware must be built, but also software written.

Some CS programs do touch DSP at the software level. An applied math program would touch even more the algorithmic side, since pretty much all DSP techniques are based on certain mathematical areas such as Fourier analysis, Complex analysis and Numerical analysis. An EE program might have a specialisation on DSP, covering both the hardware and software sides but at a more "applied", less "proof mathematical" level, than what mathematics programs do.

The field is that wide though that it might also be better to focus on only one thing and do it well, rather than think in a "jack of all trades" way, i.e. too broadly. That is, if your interest is in software and mathematics, then focus on those. If your interest is in hardware, then focus on that. Applied math is sort of a middle ground. It's both theoretical, but it will lack practice such as electronics design, even if it'd give you conceptual tools for understanding electromagnetism and digital logic.


Fields of math needed for digital filter design I want learn digital filter design. My knowledge of math is at high school level. I can learn math through the Internet. Then, what fields of math do I have to learn?

To get started:

Complex numbers

The frequency response of a filter is easier to understand complex-valued, describing both the magnitude frequency response and the phase frequency response. You will be able to understand poles and zeros, which can be complex. Complex numbers enable you to have negative frequencies, which will make math simpler.

Trigonometry

$\sin$, $\cos$ and their relation to the complex exponential $e^{i\alpha} = \cos(\alpha) + i\sin(\alpha)$ are important. Sinusoidal functions will be passed through filters with only their amplitude and phase affected.

Differentiation

To find at what frequency a simple filter peaks or dips, you can solve at what frequency the derivative of its magnitude frequency response is zero.

Integration

Integration is needed for the Fourier transform and the inverse Fourier transform.

Fourier transform

Fourier transform enables you to go from an impulse response to a frequency response and back. Also things you do in the time domain often have a simple counterpart in the frequency domain, and vice versa.


Two-Box-Model of a nonlinear amplifier A nonlinearity with memory can be modelled by a two-box-model, which consists of a filter and a memoryless nonlinearity. I am referring to chapter 5.3.2 of the book "Simulation of Communication Systems" by M.C.Jeruchim et. al. here's a link to Google books

I am given the AM-AM- and the AM-PM-characteristic of an amplifier. I want to model the behavior of the amplifier to a two-tone input signal. To do so, I want to set up a two-box-model. The memoryless nonlinearity should be described by the given characteristics.

But what does a filter before the memoryless nonlinearity do?
What does an ARMA filter do?
Can an ARMA filter be used to filter an analog signal (sine-waves) ?

The filter is a linear time-invariant (LTI) system that adds memory to the model. Combined with a memoryless nonlinearity, the filter allows to model non-linear systems with memory, unlike a pure AM-AM and AM-PM model, which models a memoryless nonlinearity. An ARMA-filter is the same as an infinite-impulse response (IIR) filter, which computes its output as a linear combination of weighted and delayed versions of the input and output signal. IIR filters need to be implemented by recursive structures.

An ARMA (or IIR) filter is a discrete-time (linear time-invariant) filter whose response $y[n]$ to an input sequence $x[n]$ is computed as follows:

$$y[n]=\sum_{m=0}^Mb_mx[n-m]+\sum_{l=1}^Na_ly[n-l]\tag{1}$$

where $b_m$ and $a_l$ are the filter coefficients, which are chosen such that a certain desired filter characteristic (e.g., low pass, high pass, etc.) is achieved.


How to implement a moving average in C without a buffer? Is it possible to implement a moving average in C without the need for a window of samples?

I've found that I can optimize a bit, by choosing a window size that's a power of two to allow for bit-shifting instead of dividing, but not needing a buffer would be nice.  Is there a way to express a new moving average result only as a function of the old result and the new sample?

Define an example moving average, across a window of 4 samples to be:

ma <= (a + b + c + d) / 4

Add new sample e:

ma_new <= (a + b + c + d) / 4 - (a / 4) + (e / 4)
ma_new(ma, oldest_sample, new_sample) <= ma - (a / 4) + (e / 4)

A moving average can be implemented recursively, but for an exact computation of the moving average you have to remember the oldest input sample in the sum (i.e. the a in your example). For a length $N$ moving average you compute:

$$y[n]=\frac{1}{N}\sum_{k=n-N+1}^nx[k]\tag{1}$$

where $y[n]$ is the output signal and $x[n]$ is the input signal. Eq. (1) can be written recursively as

$$y[n]=y[n-1]+\frac{1}{N}(x[n]-x[n-N])\tag{2}$$

So you always need to remember the sample $x[n-N]$ in order to compute (2).

As pointed out by Conrad Turner, you can use an (infinitely long) exponential window instead, which allows you to compute the output only from the past output and the current input:

$$y[n]=\alpha x[n]+(1-\alpha)y[n-1]\tag{3}$$

but this is not a standard (unweighted) moving average but an exponentially weighted moving average, where samples further in the past get a smaller weight, but (at least in theory) you never forget anything (the weights just get smaller and smaller for samples far in the past).


Detecting outliers/noise from sensor data I am trying to detect outliers/noise as indicated on the diagram below from sensor data. Can anyone advice how to go about it? I can only do this in python, so are there libraries in python that I can leverage? Is there an example that can be given. Sample data am using has timestamps and the value.
Also how can I also get the frequency of the noise/outlier occurence?

Regarding your next question, I will try to give you some long description explainig a general way of transforming filters defined by their transfer function into discrete domain. If you are not interested in the details, you can skip to the bottom, where I wrote a C code (tested).

The transfer function of a first-order low pass filter (with unit gain) is defined as follows:
$$K(s) = \frac{1}{1+sT}$$
And its angular cutoff frequency defined as $$\omega_c = \frac{1}{T}$$
Therefore the cutoff frequency expressed in Hertz becomes
$$f_c=\frac{1}{2 \pi T}$$
$$T=\frac{1}{2 \pi f_c}$$
Whatever goes above this frequency, will be cut off with gain equal to -3dB/decade for the first order low pass filter.

Either using bilinear transform or by expressing the gain in form of differential equations and then discretizing them with Euler scheme you can obtain a discrete domain equation that expresses the gain. I prefer bilinear transform as it is more straightforward.
In order to apply bilinear transform we perform the following substitution:
$$ s = \frac{1}{T_s} \frac{1-z^{-1}}{1+z^{-1}}$$
Where Ts (not to be mistaken with T) is sampling period, so in your case it is the number of seconds between each sample. I.e. for 44100Hz sampling frequency, Ts becomes 1/44100 seconds.

Now, after performing this substitution and simplifying terms, the discrete time transfer function (Z transform) of your filter becomes:
$$\frac{1+z^{-1}}{(1+2T_s/T) + (1-2T_s/T)z^{-1}}$$
Now lets recall the definition of Z transform. In general form it is expressed as the following:
$$ H(z) = \frac{b_0 + b_1z^{-1}+ b_2z^{-2} + ... + b_Mz^{-M}}{1 + a_1z^{-1}+ a_2z^{-2} + ... + a_Nz^{-N}} = \frac{Y(z)}{X(z)}$$
Where bare numerator coefficients from 1 to M (highest power) and a are denominator coefficients.

Now, a difference equation can be derived in the following way:
$$y(n) = b_0x(n) + b_1x(n-1) + b_2x(n-2) + ... + b_Mx(n-M)  - a_1y(n-1) - a_2y(n-2) - ... - a_Ny(n-N) $$
Where n means n-th sample, y is your output and x is your input. After applying the above, the simplified form of difference equation for your filter becomes:
$$y(n) = (1-\alpha)y(n-1) + \alpha x(n)$$
Where
$$ \alpha = \frac{T_s}{T + T_s}$$

This means that n-th sample (the current one) in filtered signal, called y(n) here, equals one minus alpha times the previous sample in filtered signal plus alpha times n-th sample in your source (noised) signal. In C language:

double lowpass (double x, double y_prev, double alpha)
{
    return (1.0-alpha)*y_prev + alpha*x;
}

int main()
{
    double x[10] = {1,2,3,4,5,6,7,8,9,10}; //signal with noise
    double y[10]; //filtered signal (will be generated)
    double fc = 500.0; //500Hz cutoff frequency
    double Ts = 1.0/44100.0; //for example
    double T  = 1.0/(2.0*3.1415*fc);
    double alpha = Ts/(T+Ts);

    y[0] = x[0]; //initialize first y sample as equal to first x sample

    printf("Input\tOutput\n");
    int i;
    for(i=1; i<10; ++i)
    {
        y[i] = lowpass(x[i], y[i-1], alpha);
        printf("%.4f\t%.4f\n", x[i], y[i]);
    }
}

Without complicating this code, for higher order filtering, you can swap y with x once it is filtered and repeat what's in the loop. It's slow and naive, but otherwise, you'd have to take a power of the transfer function (i.e. K(s)^2 for second order filter) and derive the discrete equation again. Thankfully with some practice, you will see a pattern in this and (even better!) you can find ready derivations of higher order low pass filters in the Internet :)

Sources:
https://en.wikipedia.org/wiki/Low-pass_filter
https://en.wikipedia.org/wiki/Bilinear_transform
https://ccrma.stanford.edu/~jos/fp/Z_Transform_Difference_Equations.html


Inverse DFT: Is there a valid / intuitive interpretation of results for non-integral timestamps? I implemented the plain DTF / inverse DFT algorithm in C++ in order to help me understand the method. As a sample input I considered
f(x) = sin( x * PI / 5 )
and collected samples for x = 0, ..., 9. So my input data was (0, 0.58778, 0.951057, ...).
Computing the DFT formula gives two non-zero coefficients: c_1 = (0 - 5i) and  c_9 = (0 + 5i).
Now, I compute the inverse DFT. For x=1, these coefficients multiplied with their basis functions give
(1/10) * exp( 2 * PI * i * 1/10 ) * (0 - 5i) = (0.293893, -0.404508i)
and
(1/10) * exp( 2 * PI * i * 9/10 ) * (0 + 5i) = (0.293893, 0.404508i)
respectively which is correct as the real parts sum up to 0.58778 (thus, reconstructing the input value) and the imaginary parts sum up to 0.  This works great for integral timestamps x=0, 1, 2, ...
My question: ... but how to interpret the results if I plug in x=0.5 (result 0.000000 - 0.951057i) or x=1.5 (result 0.000000 - 0.587785i) for example? I see that it is not DFT's promise to provide a good interpolation here. But is there an intuitive interpretation?  (Came across this because I wanted to plot the DFT/IDFT result as a time-continuous graph and don't know how to deal with the complex results between)
EDIT
I add my C++ code here:
void dft( const std::vector<double> & vInput, std::vector< std::complex<double> > & vCoeff )
{
    const std::complex<double> i( 0., 1. );

    const int nElem = vInput.size();

    vCoeff = std::vector< std::complex<double> >( nElem, std::complex<double>(0., 0.) );

    for( int k = 0; k < nElem; ++k )
    {
        for( int j = 0; j < nElem; ++j )
        {
            vCoeff[k] += exp( -2. * M_PI * i * std::complex<double>( j, 0. ) * std::complex<double>( (double)k / nElem, 0. ) ) * vInput[j];
        }
    }
}

std::complex<double> inverseDft( const std::vector< std::complex<double> > & vCoeffs, const double t )
{
    std::complex<double> result( 0., 0. );
    const std::complex<double> i( 0, 1 );
    const int nElem = vCoeffs.size();

    for( int j = 0; j < nElem; ++j )
    {
        result += vCoeffs[j] * exp( 2. * M_PI * i * std::complex<double>( j, 0. ) * std::complex<double>( t / nElem, 0. ) );
    }

    return (result / nElem);
}

This is how I define my input:
std::vector<double> vInput;

for( int i = 0; i < 10; ++i )
{
    vInput.push_back ( sin( i * M_PI / 5 ) );
}

And this is the output I did not expect for t=0.5, t=1.5, etc.:
inverseDft( vCoeffs, 0. );  // okay (0.000000, -0.000000)
inverseDft( vCoeffs, 0.5 ); // ?? (0.000000, -0.951057)
inverseDft( vCoeffs, 1. );  // okay (0.587785, -0.000000)
inverseDft( vCoeffs, 1.5 ); // ?? (0.000000, -0.587785)

EDIT2
Great to have these expert answers here!  I can contribute the 3D plot which supports my intuition and therefore might also be helpful for others. It shows f(t) computed by IDFT and how it travels over the original data points (crosses). In order to achieve this, it makes use of the whole space including the complex plane.

Andy Walls provided the answer with the direct interpretation of an N-point inverse DFT with fractional time stamps. However, there is a slighty different (or totally should I've said) interpretation of a non-integer IDFT argument, which reduces to a good interpolation formula from the inverse DFT.
Consider the notation with $N,M,k,n$ being integers.
$$x_i[n] = \frac{1}{N} \sum_{k} X[k] \cdot e^{ j \frac{2\pi}{N} k (n) } \tag{0} $$
Lets say Eq.0 is a valid inverse DFT for some N-point $X[k]$, and results in a sequence $x_i[n]$, we ask what happens if we evaluate the sum at $n = n/M$; i.e., can we also plug $n/M$ into RHS $x_i[n]$ as given by Eq.1 below ?
$$x_i[n/M] = \frac{1}{N} \sum_{k} X[k] \cdot e^{ j \frac{2\pi}{N} k (n/M) } \tag{1} $$
Where $X[k]$ is the $N$-point DFT of some sequence $x[n]$, and $x_i[n]$ is the time sequence that result from the inverse transform evaluated at non-integer time indices. And the range of summation is deliberately left out.
With the index notation of $n/M$ what's implied is not very clear. Plugging $n/M$ into an inverse DFT is formally not allowed as the argument $n/M$ is meaningless for the time-domain sequence. Yet it's possible to interpret the notation as an indicator of an oversampling of the implied interpolation applied on the $x[n]$ sequence.
In that sense, I would re-write the inverse DFT as
$$x_i[n] = \frac{M}{L} \sum_{k} X[k] \cdot e^{ j \frac{2\pi}{N*M} k n } \tag{2}$$
Now $x_i[n]$ looks like an $L = M \times N$ point inverse DFT of the N-point forward DFT $X[k]$ of $x[n]$. But there're two alternative interpretations; the first one is explained by Andy Walls, in which the inverse DFT is applied to the zero appended DFT $X[k]$ which looks like:
$$x_i[n] = \frac{M}{L} \sum_{k=0}^{L-1} X[k] \cdot e^{ j \frac{2\pi}{L} k n } = \frac{M}{L} \sum_{k=0}^{N-1} X[k] \cdot e^{ j \frac{2\pi}{L} k n } \tag{3}$$
Using Matlab you can implement Eq.3 as:
xi = M*ifft( fft(x,N), L );

This is a terrible interpolation, as explained in the other answer.
But I would call it terrible because of the miss-placed zeros in the inverse DFT formula. It performs zero padding by appending zeros to the tail of $X[N-1]$ until the end $k=L-1$
Now, the other alternative is to add the zeros into the middle of $X[k]$, with the following associated inverse DFT, for $N$ even :
$$x_i[n] = \frac{M}{L} \big( \sum_{k=0}^{N/2} X[k] \cdot e^{ j \frac{2\pi}{L} k n } +\sum_{k=N/2+1}^{L-N/2} 0 \cdot e^{ j \frac{2\pi}{L} k n } +  \sum_{k=L+1-N/2}^{L-1} X[k] \cdot e^{ j \frac{2\pi}{L} k n } \big) \tag{4}$$
In Eq.4, the orignal $N$-point DFT $X[k]$ is repeated $M$ times, (indicator of expanding $x[n]$ by $M$), and windowed by an $N$ point rectangular window with a symmetric placement (representing the frequency domain effect of an ideal lowpass brickwall filter applied on the expanded time-domain sequence), and then an L-point inverse DFT is performed to get the $L$-point interpolated sequence $x_i[n]$. This can be achieved in MATLAB with the following line :
xif = M*real(ifft([X(1:N/2+1), zeros(1,L-N),X(N/2+2:N)],L));

Now this is a quite good interpolator of $x[n]$ at a rate of $M$ times oversampling it...


Time-averaging as a low-pass filter? I run into the following filter that is supposed to do some type of time averaging of $x$

$$y[n]=G\cdot y[n-1]+(1-G)\cdot\dfrac{x[n]-x[n-1]}{T} $$

where

$$G=e^{-2\pi c T} $$

$T$ is a time interval, and the samples of $x$ and $y$ come every $T$ seconds (we can assume $y[0] = 0$).

The definition includes a quantity $c$ that I think acts as some sort of cut-off frequency, but I am not sure.  Can this type of time-averaging be seen as a low-pass filter? If so what type of filter, and how can I relate $c$ to $T$ considering that the samples come every $T$ seconds?

first of all, if i were the OP, i would fix the notation a bit.  $x(t)$, $y(t)$ are continuous-time signals and $t\in\mathbb{R}$ is a continuous argument.  $x[n] \triangleq x(nT)$ and $y[n] \triangleq y(nT)$ are discrete=time signals and and $n\in\mathbb{Z}$ is a discrete argument, an integer.  square brackets mean that an integer goes into them.  the normal convention is that $T=\frac{1}{f_s}$ is your sampling period.

a simple moving-average filter with equal-weighting in the average is

$$ \begin{align}
   y[n] & = \frac{1}{N} \sum\limits_{i=0}^{N-1} x[n-i] \\
        & = \frac{1}{N} \left(x[0] + \sum\limits_{i=1}^{N} x[n-i] - x[n-N]\right) \\
        & =  \frac{1}{N} \sum\limits_{i=1}^{N} x[n-i] + \frac{1}{N} \left(x[0]  - x[n-N]\right) \\
        & =  \frac{1}{N} \sum\limits_{i=0}^{N-1} x[n-1-i] + \frac{1}{N} \left(x[0]  - x[n-N]\right) \\
        & =  y[n-1] + \frac{1}{N} \left(x[0]  - x[n-N]\right) \\
\end{align} $$

that is representing both the non-recursive implementation and the recursive.  the recursive implementation works as long as you subtract exactly what you earlier added, something that might not happen in a floating-point context.  both are FIR, even though the recursive implementation is something we would call a Truncated IIR,  still an FIR.

to get the frequency response, apply the Z-Transform to either form, and solve for $Y(z)$ in terms of $X(z)$.  then compute the transfer function:

$$ H(z) \triangleq \frac{Y(z)}{X(z)} $$

then substitute $z = e^{j \omega}$  and solve for the $\omega$ such that $\left|H(e^{j \omega})\right|^2 = \frac{1}{2}$.  that $\omega= 2\pi f T$ is your -3.01 dB cut-off frequency.

try doing that math yourself.  we'll be here to help guide you.


What mother wavelet should be chosen as an alternative to the STFT? The Short Time Fourier Transform (STFT) is used to identify time localized frequency content of a signal. The STFT operates by chunking an input signal into blocks and performing FFT on the block, and then recording the frequency coefficients for that point in time. However, the STFT enforces fixed time and frequency resolution which is not always desired.
A wavelet transform is an alternative method for identifying the frequency content of a signal where the window size typically varies based on the frequency, thus creating an inverse relation ship between time localization and frequency localization. That is: high-frequencies have good time localization but poor frequency localization, and low-frequencies have good frequency localization but bad time localization.
Let's say I want to use a wavelet transform as an alternative to the STFT. What specific wavelet should I use? Do some wavelets work better for encoding frequency content?

Just because the math is easier, I might recommend the Morlet or possibly the Chirplet as the mother wavelet.


Radar signal processing flow chart I'm new to Radar in general, and I'm studying RADAR using online resources.  I came across the following slide which shows the process flow of the Radar signal processing.

The thing I don't get in the figure is the signal processing block and the general purpose Computer block. The following is my basic understanding of the process flow which is very different than the flow in the picture:

Pulse compression -> matched filter -> Threshold -> Range Detection -> Clutter rejection (Doppler processing) -> Doppler Detection.

The reasoning is that we need to find the target's range gate before we can process the target's Doppler information, and to detect the target's range, we need to set some threshold to separate our target from other clutters and noise. After we locked the target into certain range gate, we can process the target's Doppler information from the I-Q data obtained after processing the received echo. Any clutters in the same range gate can be rejected by using a set of Doppler filter banks as the target and clutters have their respective Doppler characteristics/frequencies.
Could someone explain to me where was my misunderstanding? I just got started, so there must be many things I got misunderstood. So thank you for your time, patience, and help.

regards,
zhengguo Sun

Part of your misunderstanding comes from the fact that there are many ways in which the radar signal processing chain is implemented. Depending on the type of radar, targets of interest, hardware, etc., some methods are more appropriate than others. We will consider pulsed-Doppler radar here.

In the chain you describe:

In modern pulse-Doppler systems using pulse compression, the act of performing pulse compression is virtually synonymous with matched-filtering. You can combine these two terms to mean the same thing and simplify your thinking. You are requiring that range detection happen before processing Doppler, but this not be the case. You can very well detect a target using its Doppler information only while being denied range. It is mainly during acquisition of the target that being denied either of the measurement dimensions is truly detrimental. In the tracking phase, you can "track through" and continue to update target velocity information using Doppler only to then later be able to receive good range measurements.

In addition, the system you describe is a more classical approach of using a Doppler filter bank to separate clutter and to determine the target's radial velocity. Many radar systems require that a signal make it past the Doppler filter bank to even consider performing a measurement. Hence why the block diagram you show has the "detection" block coming after the "Doppler filtering" block.

What the block diagram suggests:

The block diagram can represent two main schemes:

Single echoes from targets that make it past the filter bank are considered to perform a measurement. If a target appears at the output of one of these filters (i.e. adequate SNR and mitigated clutter), you can declare a detection where range and Doppler are measured simultaneously. Multiple pulses may be integrated to improve SNR.
The radar system performs Doppler processing by forming a range-Doppler map. This map can be seen as a 2D array formed by collecting multiple target returns after matched filtering for multiple pulse-repetition intervals (PRIs). Once a certain number of matched-filtered returns are collected (usually a power of 2), a DFT is performed in what is called the "slow-time" dimension. In other words, a DFT is performed across every range gate (a range gate is equivalent to each sample of the matched-filter output). After this step, you now have a 2D map where one axis is the delay (via the matched-filter) and the other axis is the Doppler frequency. Detection now finds the cell(s) according to a certain threshold. A cell is "detected" if it passes a certain threshold, and you have performed range and Doppler measurements simultaneously by virtue of simply picking the cell. Hence, you use the range-Doppler map to discriminate using Doppler and perform detection at the same time! You can go further and use both measurement dimensions to discriminate possible detections. One application being to keep track of an airborne radar's altitude return.

These are just the basics! You can implement all sorts of filtering and signal processing algorithms on the range-Doppler map alone to get even more information about your environment and your targets of interest.


How to soft decode DQPSK? I'm successfully soft-decoding D-BPSK by taking the dot-product of the constellation-position of the symbol and of the previous symbol.  If the result is >= 1, then the symbol phase hasn't changed and the bit is a zero.  If the result is <= -1 then the phase has shifted and the result is a one.  In between -1 and 1 the result is a soft 0 or soft 1.

I can't figure out how to do the same thing with D-QPSK.  I can use just the phase, but this throws away a lot of information that could help the soft-decoder.

This paper explains how to do it and gives a formula (10):

$b_1 = \mathrm{Re}\{s_n s^*_{n-1}\}, b_2 = \mathrm{Im}\{s_n s^*_{n-1}\}$

But I don't understand the notation — what does a * floating above mean?  I tried just multiplying the complex numbers and taking the real and imaginary parts but this didn't work.

Since the constellation can rotate, how can the two axes be teased apart?

Two successive symbols in the demodulator are $Z_1 = (X_1,Y_1)$ and $Z_2 =(X_2,Y_2)$
where $X$ is
the output of the I branch and $Y$ the output of the Q branch of the receiver. The
hard-decision DBPSK decision device considers the question:

  Is the new symbol $Z_2$ closer to the old symbol $Z_1$ or to the
  negative $-Z_1$ of the old symbol?

and thus compares

$$(X_2-X_1)^2 + (Y_2-Y_1)^2 \gtrless (X_2+X_1)^2 + (Y_2+Y_1)^2$$

which can be simplified to a sign comparison on $\langle Z_1,Z_2\rangle = X_1X_2+Y_1Y_2$.
Note that this is essentially asking

  Are the two vectors $Z_1$ and $Z_2$
  are pointing in roughly the same direction (in which case the inner product
  or dot product is positive) or in roughly opposite direction (in which case
  the dot product is negative)?

A third viewpoint thinks of $Z_1$ and $Z_2$ as
complex numbers and asks

  Is $\text{Re}(Z_1Z_2^*) = X_1X_2+Y_1Y_2$ positive or negative?

The soft decision decision device simply
passes on the exact value of the dot product to the soft decision decoder
which may opt to quantize dot products that are very large in magnitude
into hard decisions and continue waffling on the rest.  This is what
the decision rule stated in the OP's question is, where large is taken as
exceeding $1$ in magnitude.

In DQPSK, the encoding uses one of two conventions:

the signal phase is delayed by $0, \pi/2, \pi, 3\pi/2$ according
as the dibit to be transmitted is $00, 01, 11, 10$
the signal phase is advanced by $0, \pi/2, \pi, 3\pi/2$ according
as the dibit to be transmitted is $00, 01, 11, 10$

Note that a DQPSK signal is not the sum of two DBPSK signals
modulated on phase-orthogonal carriers, but the I and Q bits
jointly affect the net carrier phase.

For demodulating a DQPSK signal, the decision device needs to ask

  Which of the four symbols $Z_1,\quad jZ_1 = (-Y_1,X_1),\quad -Z_1,\quad -jZ_1 = (Y_1,-X_1)$
  is $Z_2$ closest to?

Thus, in addition to the comparison

$$(X_2-X_1)^2 + (Y_2-Y_1)^2 \gtrless (X_2+X_1)^2 + (Y_2+Y_1)^2$$

it is necessary to compare

$$(X_2+Y_1)^2 + (Y_2-X_1)^2 \gtrless (X_2-Y_1)^2 + (Y_2+X_1)^2$$

which works out to looking at $\text{Im}(Z_1Z_2^*)$ in
addition to $\text{Re}(Z_1Z_2^*)$ and making the decision
according as which quantity has the largest magnitude
and the sign of the largest magnitude.  The details
of how the soft-decision decoder uses the decision
statistic $Z_1Z_2^* = (\text{Re}(Z_1Z_2^*), \text{Im}(Z_1Z_2^*))$
will determine how these numbers are further
massaged.


Butterworth filter transfer function in time domain Preliminaries
I need to know the response of Butterworth filter in time domain of the following input (at least numerically):
$$
x(t) = \theta(t)(1-\exp\{-at\})\exp\{-bt\}, \quad a,b>0,\quad a\gg b
$$
In order to find the response I decided to find the transfer of Butterworth filter in time domain and then convolve it with $x(t)$ numerically.
But the result I've got doesn't seem to be correct. Here how it went.
Transfer function in $s$-domain
Transfer function of Butterworth filter is the following:
$$
H(s) = \prod_{k=1}^{n}\frac{1}{(s-s_{k})}, \quad s_{k} = \omega_{c}e^{i\frac{2k+n-1}{2n}\pi}
$$
So in order to find transfer function in time domain I used this equality:
$$
h(t) = \sum_{s_{k}}\mathrm{res}H(s_k)\exp\{s_k t\}
$$
Transfer function in $t$-domain
First of all, I found that $\Im{h(t)} = 0$ i.e. it is real. I used this fact as a cross check. I wrote the code that calculates $h(t)$ using complex numbers. I found the following result:

The question is : could somebody verify if it looks right? I haven't found any pictures of the Butterworth transfer function in time domain.
P.S
I could post the code I used to produce those results, but keep in mind it is based on ROOT CERN. You can see the code and tests here.
Update
I have performed one more test: calculation of the magnitude of the frequency response of Butterworth filter of the order $n=3$ and cut frequency $\omega_{c}=2$ rad/s. The result is perfect:

A plot of the normalized impulse responses, for the n = 2 through 10 Butterworth low pass filters, are given by H.J. Blinchikoff, A.I. Zverev, "Filtering in the Time and Frequency Domains", Wiley-Interscience, John Wiley & Sons, NY, ©1976, p. 113. This is shown below. They do not give the h(t) expressions in the book, at least where I have looked thus far. They say the figure is from Zverev's earlier book: A.I. Zverev, "Handbook of Filter Synthesis", Wiley, NY, ©1967.

EDIT: Looking back in my old computation notebooks, I found the impulse expression for the 3rd order Butterworth low pass filter. This pic is from 1983:

Hopefully, this is algebraically the same as your result.

EDIT 2: Gregg 1 gives $h(t)$ for the $M = 2$ Butterworth LPF:
$$h(t) = e^{-x}sin(x)U(t)$$
where $x = t/τ√2$, $τ = RC = 1/(2πB)$ and $B$ is the corner frequency in Hz. Butterworth LPFs are fine for many purposes, but not particularly good for others, e.g., they have relatively poor 1% settling times with step inputs 2.

References:

W.D. Gregg, “Analog and Digital Communication”, Wiley, NY, ©1977.
E. Voigtman, J.D. Winefordner, “Low-pass filters for signal averaging”, Rev. Sci. Instrum. 57 (1986) 957-966.

EDIT 3: Even though Matt has provided an excellent answer, for completeness I decided to test my 1983 $h(t)$ expression (see above) by comparing it with the numerically generated impulse response from a computer simulation $H(s)$ transfer function block. As per the OP, I assumed $2πB = 1$. The next figure shows that the results are the same.


What does opencv `mean` function do when passing a mask of zeros? INTRO

I am looking at this project to implement some changes.

The project is written in C++ using openCV library. Its goal is to detect cars in parking area to calculate the number of free spots. It does this in real-time (takes the input of a camera and processes it frame by frame).

The detection is performed on areas of the image that are rectangles whose vertices are loaded from a file. For each park (i.e. rectangle) to be monitored there is an instance of the class Parking defined in the files Parking.*** and all of them are collected in a vector (named parking_data in the following lines of code)

CODE

As far as I understood the processing starts on line 83 of the file main.cpp. The interesting lines are these:

cv::cvtColor(frame, frame_gray, cv::COLOR_BGR2GRAY);
cv::GaussianBlur(frame_gray, frame_blur, blur_kernel, 3, 3);
//...
for (Parking& park : parking_data)
{
  // Check if parking is occupied
  roi = frame_blur(park.getBoundingRect());
  cv::Laplacian(roi, laplacian, CV_64F);
  delta = cv::mean(cv::abs(laplacian), park.getMask());
  park.setStatus( delta[0] < atof(ConfigLoad::options["PARK_LAPLACIAN_TH"].c_str()) );

cvtColor converts the frame to grayscale then a gaussian blurring filter is applied. Then it iterates to the elements of the vector parking_data containing the instances of the class Parking which represent the area in the image corresponding to the parking spots.

getBoundingRect should return the vertices of the rectangle so to get a subimage of the blurred frame. Then a Laplacian filter is applied.

QUESTION

The line I do not understand is delta = cv::mean(cv::abs(laplacian), park.getMask());.

I understand it is calcualating a mean of the absolute value of the laplacian subimage. The mask is defined only in line 49 and used in the followings of the file Parking.cpp as

mask = cv::Mat::zeros(bounding_rect.size(), CV_8UC1);
for (cv::Point p : contours_points.at(0))
{
  polygon_points_in_bounding_rect.push_back(cv::Point(p.x - bounding_rect.x, p.y - bounding_rect.y));
}
vector<vector<cv::Point>> contours;
contours.push_back(polygon_points_in_bounding_rect);
cv::drawContours(mask, contours, -1, cv::Scalar(255), CV_FILLED);

so it is a matrix of 0s.
Then the value is compared to a given threshold.

OpenCV on the function says about the mean function.

When all the mask elements are 0's, the functions return Scalar::all(0)

So what is the mean function calculating when passing a mask of zeros?

The function cv:mean calculates the mean value of the marked out area in the InputArray as you said. So the function looks at the area which is filled with the mask (zeros) and calculates the mean value of the InputArray which is not filled with zeros.

I know the code you are showing here. I don't understand why the function cv:Laplacian is used here because in my opinion this is not relevant in this case. Maybe someone else knows why he used the Laplacian array instead of using the blured_gray array.

I hope I could help you.


Do all Harmonics of a periodic sound start from a 0 base I've been thinking about additive synthesis and I'm wondering if all sounds begin at (theoretically) phase 0?

It's common to not start all harmonically related sinusoids with a phase of 0 at some initial time point.  Varying the phase relationships helps reduce the peak-to-average power ratio, which allows a higher average volume without clipping or otherwise exceeding some peak limit.


Energy and SpectralFlatnessMeasure taken for PCM audio VoiceActitivyDetection I'm trying to implement the VAD algorithm in this paper:

http://ms12.voip.edu.tw/~paul/Papper/Steganography/iLBC/(VAD)Real-Time_VAD_Algorithm.pdf

I'm polling a microphone on an Android phone, and getting short[]'s which I want to perform operations on.  The audio data is 16bit PCM data.

I've gotten the basics of the algorithm done, but there are two issues I don't really understand yet.

1)  Step 3-1- of the above algorithm says to

  "compute frame energy".

I'm not entirely sure what this means.  Should I take the PCM data and convert it to Db's? Should I generate this off of the magnitude values derived from an FFT?

2) In section 2 (1) of the paper, There is an equation for Spectral Flatness given:

SPM = 10log10(Gm/Am)  or Spectral Flatness Measure = 10*Log10(GeometricMean/ArithmeticMean)

There is a brief description underneath which reads:

  "Where Am and Gm are the geometric means of speech spectrum
  respectively".

I understand what the general terms of this equation are, however I don't understand what I should be summing to generate the Gm and Am terms.  At first I thought it was the energy values computed in part 1) of my question, however rereading it seems to instead relate to the energy present only in the speech spectrum portion of a given frame.  If thats the case what range should I use and how do I generate the value to sum/multiply?

The frame energy can be computed in the time domain as well as in the frequency domain. For a frame of $N$ samples, due to Parseval's theorem, you get for the frame energy

$$\sum_{n=0}^{N-1}|x[n]|^2=\frac{1}{N}\sum_{k=0}^{N-1}|X[k]|^2$$

where $x[n]$ are the time domain samples, and $X[k]$ are the complex DFT coefficients.

The spectral flatness is commonly computed from the power spectrum, or from the squared magnitudes of the DFT coefficients. I can't open the link to the paper, but unless otherwise stated, I would compute the means using all values $|X[k]|^2$ in the current frame.


Fitting piecewise splines to noisy data I have a system that gives me a noisy data set similar to the one generated by this matlab/octave code. The y-axis represents the signal intensity and the x-axis represents spatial distance.

noise = randn(1,512);
points = gaussian(512,.01);
noisy = points + .02*(noise');
plot(noisy)

The actual signal will be asymmetric and not the result of a gaussian function(but close). I can't create the noise accurately by hand and the underlying signal algorithmically so image these two together.

The goal is to construct an interpolating curve, using bicubic splines, based on only one data set of this signal. The signal cannot be captured multiple times and averaged together to remove the noise. Since the signal is noisy, I can't simply pick points every N samples as the interpolation points so I need to do some signal conditioning first.

I will not be able to use an exact model of the underlying physical process but the sources of noise are known and I can make some assumptions regarding the maximum slope and SNR. The interpolation is used to correct for manufacturing limitations so the curve is intended to be smooth(exclusive of noise) but for cost reasons is not. Therefore we don't expect any large deviations from this example.

My initial approaches used an M-point moving average filter to smooth the signal, then used polyfit handle the ends of the signal. This worked but could potentially reduce the center peak, if M is large enough. I've also since learned smoothing signals before interpolation is generally discouraged.

What would be a better approach to achieve this?

It seems like you don't want to make any prior assumptions about the exact shape of the signal, other than assuming that it is ``smooth''. This is best attacked using a non-parametric function fitting/regression approach. Use either
a) spline fitting or b) kernel smoothing.

Two reasons why this is a good idea:

These methods are optimal in the sense that no other non-parametric filtering technique (moving average, etc.) can give you a better worst case error decay rate. This answers your main question: there is no better method if you care about minimax error.
Spline fitting and kernel regression codes are very easily available.

A quick Google search gives me the following:
In Octave interp1 has a 'cubic' option:
http://www.gnu.org/software/octave/doc/interpreter/One_002ddimensional-Interpolation.html

Open source code for kernel smoothing:
http://code.google.com/p/pmtk3/source/browse/trunk/matlabTools/stats/kernelSmoothingRegression.m?r=2645&spec=svn2776


Sparse Inverse Fourier Transforms I'm looking to compute an IDFT but my input is very sparse:
Example:

IDFT length: 7,408,800 (complex floats)
Sparsity: 96.61% to 99.99%

I found this Sparse FFT website but it looks like the library deals with sparse inputs to standard DFT - not IDFT - also it doesn't look like it is not for consumer use.
Are there any techniques I can employ to improve performance for sparse IDFT - before even changing my DFT library? And now, are there any available libraries that give performance improvements for spare IDFT?

The FFT is faster than the naive DFT through matrix-vector products because it can reuse many intermediate results.
However, for inputs this sparse, there's really not much that can be re-used, even in the best case.
So, the most efficient way here is probably to work straight forward:

Allocate the space ($N=7408800$ values) for your output (initialize with zeros) $o$.
Go through your input $i$. Nested loops:

For every nonzero input index $n$:

For every output index $k$: add $i_ne^{j2\pi \frac{nk}N}$ to $o_k$.

Notes:

since with sparse input, your output is non-sparse, it's usually faster to traverse the input only once, but the output many times, instead of the other way around (that is because memory is much faster when addressed linearly)
to calculate the complex exponential, you can use SIMD-accelerated calculations, because the values are in consecutive addresses and have proportional phase. LibVOLK has a number of implementations that might be helpful for fast calculation here.


MSK spectral efficiency explanation Why is MSK (minimum shift keying or FSK with a modulation index of 0.5) said to be a spectrally efficient form of FSK?

In traditional FSK, the frequency separation between pulses is $1/T$, where $T$ is the symbol rate. This guarantees orthogonal pulses while maintaining a continuous phase. The phase is continuous because a pulse with the next higher frequency traverses exactly one more cycle than the one with the lower frequency (during one symbol interval).

However, you can easily show that in order to guarantee orthogonality between pulses, a frequency separation of $1/2T$ is actually sufficient. This is indeed the minimum frequency separation guaranteeing orthogonality. Now a higher frequency pulse traverses only half a cycle more than the pulse with the next lower frequency. The smaller frequency separation reduces the bandwidth compared to traditional FSK, whereas the probability of error remains unchanged due to the orthogonality of the pulses. This is why MSK is often referred to as bandwidth-efficient (compared to traditional FSK).

Note that for MSK to have a continuous phase the initial phase of a pulse must be either $0$ or $\pi$, depending on the phase of the previous pulse (because some pulses have a phase value of $\pi$ at the end of a symbol interval due to the minimum frequency separation of half a cycle per symbol interval).

This additional "phase encoding" inherent in MSK can actually be used in the receiver by applying the Viterbi algorithm. In this way, MSK can perform better than traditional FSK. It can overcome the 3dB loss of FSK as compared to PSK.

Also note that MSK can be viewed and implemented as a special case of offset keyed PSK (with sinusoidal pulse shapes).